Introducción
El análisis de series climáticas es siempre un desafío y más aún en regiones donde la disponibilidad de las mismas es escasa. Este déficit es un problema importante para la realización de análisis de riesgo o de modelación que demandan una gran cantidad de datos. Los generadores estocásticos de clima producen series diarias de variables climáticas con propiedades estadísticas consistentes a los de los datos observados. Estas series sintéticas son un insumo muy importante para el análisis cuantitativo del riesgo de desastres naturales. Este tipo de aplicaciones requiere de largas series temporales para poder capturar la variabilidad natural del clima de una región. Las principales variables de importancia son temperaturas máxima y mínima y precipitación. En muchas regiones estos datos se encuentran incompletos o son inexistentes. Los registros suelen tener una duración insuficiente o sólo estar disponibles datos agregados mensualmente. Para suplir este déficit, los generadores estocásticos pueden generar series que satisfagan las necesidades.
La mayoría de los enfoques tradicionales para la generación de series estocásticas están limitados por su capacidad de generar datos solamente en localidades para las que se cuenta con observaciones (por ejemplo, donde existen estaciones meteorológicas). Otra desventaja de algunas de estas herramientas (sobre todo los generadores no paramétricos basados en remuestreo de observaciones) es que solamente pueden producir valores dentro del rango observado en el registro histórico. En este proyecto, el generador desarrollado se basó en el modelo diseñado por Verdin et al. (2016). Este generador estocástico diario multivariado produce series sintéticas de precipitación y temperaturas máxima y mínima. El generador de Verdin et al. (2016) utiliza cuatro modelos estadísticos para modelar la ocurrencia y montos de precipitación y temperatura máxima y mínima basados en Modelos Lineales Generalizados (GLMs). Estos modelos son paramétricos y utilizan regresiones lineales entre las variables para modelarlas. El proceso de ocurrencia de precipitación se modela como una regresión probit mientras que los montos se ajustan a una distribución aleatoria Gamma. Las temperaturas máximas y mínimas son modeladas como variables autorregresivas condicionadas por la precipitación. El generador produce series espacialmente correlacionadas tanto para estaciones individuales como en grillas regulares. Una de las limitaciones fundamentales del generador de Verdin et al. (2016) es que las temperaturas máxima y mínima se generan en forma independiente para cada día, por lo cual la amplitud térmica diaria simulada a veces no es realista. Los diagnósticos previos realizados en base al generador de Verdin et al. (2016) demostraron que algunas propiedades de las series sintéticas producidas no reflejaban fielmente características importantes para el análisis de las sequías, por ejemplo, la persistencia de secuencias de días secos y lluviosos. Por este motivo, en este trabajo solo se mantuvo la estructura general del modelo de Verdin et al. (2016) y se modificaron todos los algoritmos (modelos estadísticos) para obtener series sintéticas más consistentes con los registros históricos y para optimizar el proceso de cálculo.
Fundamento
A partir de los datos climáticos históricos, comienza el proceso de ajuste de un modelo estadístico que permita representar el comportamiento de cada una de las variables para cada localidad. El generador está dividido en cuatro modelos aditivos generalizados: dos para modelar la precipitación y dos para las temperaturas máxima y mínima, respectivamente. El concepto principal detrás de la generación de series sintéticas es que cada valor de una variable puede ser considerado como la suma de una componente climática (clima local) y otra componente meteorológica aleatoria o tiempo local (Kleiber et al., 2013), es decir
\[ X_{i,s} = clima\ local + tiempo\ local \] donde \(X_{i,s}\) corresponde al valor de la variable X en el día i en la localidad s; clima local corresponde a un valor medio de la variable para el día i en la localidad s y tiempo local corresponde a un estado particular de la atmósfera en el día i en la localidad s. El componente climático es especificado a través del ajuste de cada uno de los cuatro modelos aditivos del generador. El componente meteorológico corresponde a los residuos de los modelos (la diferencia entre los valores históricos y el componente climático estimado), es decir, a la variabilidad no explicada por los mismos. El proceso de ajuste culmina con los parámetros ajustados para los cuatro modelos mencionados. Posteriormente, se generan datos para los años a simular que corresponden a los valores medios de la distribución para cada día del año (clima local). Luego, se simulan una serie de valores aleatorios, a partir de los residuos de cada modelo, que corresponden a la variabilidad propia de cada realización (tiempo local).
Figura 1. Componentes clima local y tiempo local
Tipos de series sintéticas
El generador GAMWGEN produce distintos tipos de series sintéticas. Desde el punto de vista temporal, las series pueden ser no condicionadas (Figura 2.a), es decir, puramente estacionarias; pseudohistóricas (Figura 2.b), que copian la variabilidad de baja frecuencia y los cambios en la serie climática observada y condicionadas (Figura 2.c) que son series forzadas a seguir la trayectoria de una salida de un modelo de cambio climático o una trayectoria arbitraria definida por el usuario. Las siguientes Figuras muestran distintos ejemplos de lo mencionado.
Figura 2.a. Tipos de series sintéticas: Series estacionarias
La línea negra corresponde a la serie temporal observada mientras que las distintas líneas celestes corresponden a realizaciones del generador. Al tratarse de series puramente estocásticas, son independientes entre sí.
Figura 2.b. Tipos de series sintéticas: Series pseudihistóricas
La línea negra corresponde a la serie temporal observada mientras que las distintas líneas celestes corresponden a realizaciones del generador. En este caso, el generador fue forzado a seguir la trayectoria en los datos observados, por lo tanto, las series generadas siguen las variaciones de los datos observados.
Figura 2.c. Tipos de series sintéticas: Series condicionadas
En esta configuración el generador sigue una trayectoria arbitraria elegida por el usuario. En este caso es una trayectoria lineal por lo que la precipitación aumenta un determinado porcentaje por año de manera lineal.
Desde el punto de vista espacial se pueden generar series en estaciones meteorológicas, en puntos que no se correspondan con estaciones o en una grilla regular. La Figura 3 muestra un ejemplo de la generación en grillas sobre el Paraguay.
Figura 3. Tipos de series sintéticas: Datos grillados
El panel de la izquierda muestra la temperatura máxima del día 1 de enero de 2000 para todo el territorio del Paraguay mientras que el mapa del panel derecho muestra la precipitación diaria para el mismo día.
Metodología
Como su nombre indica, el generador está basado en Modelos Generalizados Aditivos (GAM, por sus siglas en inglés). Como se mencionó en le Introducción, este generador está basado en uno similar desarrollado por Verdin et al. (2016). Dicho generador estocástico utilizaba modelos lineales generalizados (GLM, por sus siglas en inglés). Los GLM son modelos muy interesantes ya que son facilmente interpretables aunque carecen de la flexibilidad necesaria para capturar complejos patrones como las variaciones estacionales. Por ello, en esta versión se cambiaron los GLMs por GAMS. Estos nuevos modelos están siendo muy utilizados en diversas áreas porque heredan lo mejor de los GLM pero son mucho más flexibles. En la siguiente sección se describirán brevemente los GAMs.
Introducción a los Modelos Aditivos Generalizados
Interpretabilidad vs Complejidad
Es posible crear modelos ajustados a los datos que son lineales y relativamente fáciles de interpretar y explicar. Existe una línea recta que representa la relación entre dos variables y atraviesa la nube de puntos. Es posible graficar la relación y hacer predicciones a partir del modelo ajustado. Pero lo modelos lineales no siempre representan bien las relaciones entre variables ya que las relaciones no siempre son lineales. Las predicciones no serían buenas a partir de estos modelos.
En el otro extremo del espectro hay toda una serie de modelos de tipo “caja negra” como las redes neuronales, random forest, árboles de regresión. Estos modelos son muy buenos para predicir pero son muy díficiles de interpretar y de entender que está sucediendo en el sistema. Son muy útiles para clasificar pero no sirven para entender como una variable se relaciona con el producto del modelo.
Los GAMs proveen un interesante punto intermedio ya que se pueden ajustar relaciones complejas y no lineales e interacciones pero estos modelos son explicitos y se pueden observar las relaciones entre variables y entender porque se produce un resultado determinado.
Relaciones no lineales
En general, las relaciones entre variables de la naturaleza no son lineales y adquieren patrones muy complejos. En la Figura 4 se muestra un ejemplo de datos sintéticas para ejemplificar este concepto.
library(mgcv)
set.seed(2) ## simulate some data...
dat <- gamSim(1, n=400, dist="normal", scale=0, verbose=FALSE)
dat <- dat[,c("y", "x0", "x1", "x2", "x3")]
Figura 4. GAMs: Relaciones no lineales.
Este diagrama de dispersión muestra dos variables que claramente se encuentran relacionadas pero no de manera lineal. Si se ajusta una regresión lineal simple a estos datos se obtienen los siguientes resultados.
ggplot2::ggplot(dat, ggplot2::aes(x2, y)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(method = 'lm', se = TRUE, ggplot2::aes(colour = "Lineal")) +
ggplot2::scale_colour_manual(name = "", values = c("Steelblue")) +
ggplot2::theme_minimal()
Figura 5. GAMs: Relaciones no lineales. Ajuste lineal. Tomado de (CITAR)
El modelo no es capaz de capturar las principales características de los datos.
Modelo lineal
Gráfico de residuos
residuals <- fortify(modelo_lineal)
ggplot2::ggplot(residuals, aes(x = .fitted, y = .resid)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(se = FALSE) +
ggplot2::geom_hline(yintercept = 0, linetype = 'dashed') +
ggplot2::theme_bw() +
ggplot2::xlab('Ajustados') + ggplot2::ylab('Residuos')## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Claramente los residuos tienen una tendencia y no están distribuidos aleatoriamente.
Sin embargo, al cambiar el modelo lineal por un gam.
Figura 6. GAMs: Relaciones no lineales. Ajuste no lineal. Tomado de (CITAR)
Con un GAM se pueden ajustar los modelos a través de funciones suavizadas o splines que pueden tomar casi cualquier forma. Al utilizar splines los GAMs pueden capturar diversos tipos de relaciones no lineales y es por esto que son tan flexibles y adaptables a distintos contextos.
residuals <- data.frame(fitted = modelo_gam$fitted.values,
resid = resid(modelo_gam)) %>%
tibble::as_tibble()
ggplot2::ggplot(residuals, aes(x = fitted, y = resid)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(se = FALSE) +
ggplot2::geom_hline(yintercept = 0, linetype = 'dashed') +
ggplot2::theme_bw() +
ggplot2::xlab('Ajustados') + ggplot2::ylab('Residuos')## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Al utilizar un GAM los residuos tienen una distribución casi aleatoria sin una tendencia clara como si fue el caso del modelo lineal. Este es sólo en sencillop ejemplo del potencial que tienen los GAMs para modelar complejas relaciones. Si se desea profundizar más en el tema se recomienda revisar el librop de Simon Wood (2017)
Modelación
Como se mencionó anteriormente, el generador estocástico tiene como base cuatro modelos generalizados, dos para temperaturas máxima y mínima y dos para la precipitación que modelan la ocurrencia y los montons diarios.
Figura 6. Modelos estadísticos
Partiendo de la base de datos se extraen los datos diarios de temperaturas máxima y mínima y precipitación. Los valores de temperaturas son utilizados para ajustar un GAM que modelará la componente climática. Estos modelo puede incluir o no medias trimestrales que servirán para condicionar el modelo a seguir un determinado comportamiento. Para el caso de la precipitación el tratamiento es diferente. El fenómeno se divide en dos: ocurrencia de precipitación y montos diarios. La ocurrencia se modela con GAM para toda la serie que puede o no estar condicionado por totales trimestrales de lluvia. Los montos, en cambio, son modelados a escala mensual. Es decir, se ajustan doce modelos, uno para cada mes del año. Se utiliza esta modalidad para capturar mejor las características propias del ciclo estacional de la precipitación de cada región.
A continuación se describirán cada uno de los modelos recién mencionados.
Modelación: Clima local: Ocurrencia de lluvia
El ajuste del modelo de ocurrencia comienza con la definición de un día lluvioso. Según la Organización Meteorológica Mundial (OMM), se considera como día lluvioso a aquel día con una precipitación igual o mayor a 0.1 mm. Este valor, si bien es el sugerido por la OMM, puede ser modificado por el usuario si así lo dispone. Una vez definida la ocurrencia de lluvia se ajusta el modelo. La ocurrencia de lluvia es una variable de tipo binaria, es decir, un día puede ser lluvioso (1) o seco (0), y es muy bien representada a través de una regresión probit Kleiber et al. (2012). Este tipo de regresiones se basan en procesos latentes Gaussianos, \(W_{s,t}\), que se modelan con la siguiente relación:
\[ O_{s,t} = \Pi_{\{W_{s,t} > 0|} \] Si el proceso \(W_{s,t}\) es positivo, significa que lloverá en el día t y en la estación s y a ese día se le asignará el valor 1. Si el proceso es negativo significa que no lloverá en el día t y en la estación s y ese día se le asignará el valor 0. El uso de este tipo de procesos latentes está justificado en que, en una región, la ocurrencia de lluvia en las distintas estaciones tenderá a estar correlacionada. La función media del proceso latente Gaussiano es una regresión entre variables que se expresa de la siguiente manera:
\[ O_{s,t} = (s, O_{s,t-1}, f(doy(t)), f(ST(t)), f(lon, lat)) \] donde \(O_{s,t}\) corresponde a la ocurrencia de lluvia en sitio y día determinado; \(s\) corresponde al efecto del sitio s (ordenada al origen); \(O_{s,t-1}\) corresponde a la ocurrencia del día previo que es un término autorregresivo; \(f(doy(t))\) corresponde a una función cíclica de los días del año para considerar el efecto de la estacionalidad sobre la ocurrencia de lluvia; \(f(ST(t))\) corresponde a una función suavizada de los acumulados estacionales de precipitación (solo se utiliza si se desea condicionar el modelo) y \(f(lon, lat)\) corresponde a las coordenadas del punto s (solo se utiliza en el modelo espacial). En la práctica, esta covariable se divide en cuatro, una para cada trimestre del año, asignándole el valor de 0 para los momentos fuera del respectivo trimestre. Es importante notar que para la ocurrencia de precipitación se usa un solo término autorregresivo. Este término es muy importante para la correcta modelización de las rachas secas y lluviosas mientras que el día del año incorpora la variabilidad intra-anual.
Modelación: Clima local: Montos diarios de lluvia
El modelo de montos diarios de lluvia difiere del anterior en que no se usan todos los datos de la serie, sino que se extraen de la base de datos los montos de precipitación únicamente para los días lluviosos. La intensidad de la precipitación para una localidad s y tiempo t es modelada como una variable aleatoria Gamma cuyos parámetros de forma y escala varían en el tiempo y en el espacio (Kleiber et al. (2012)). La función Gamma ha sido ampliamente utilizada en la región para la modelación de acumulados de lluvia. El modelo de montos puede ser expresado de la siguiente manera:
\[ I_{s,t} = (s, O_{s,t-1}, f(ST(t)), f(lon, lat)) \]
donde, \(I_{s,t}\) corresponde a los montos de precipitación en un sitio y día determinado; \(O_{s,t-1}\) corresponde a la ocurrencia del día previo para considerar la autocorrelación temporal; \(f(ST(t))\) corresponde a una función suavizada de los acumulados estacionales de precipitación (solo se utiliza si se desea condicionar el modelo) y \(f(lon, lat)\) corresponde a las coordenadas del punto s (solo se utiliza en el modelo espacial). A diferencia del modelo anterior – que modela los años calendarios completos a través de la estacionalidad del día del año – la serie de montos diarios es dividida en función del mes del año para ajustar una distribución a cada mes para obtener así parámetros mensuales más precisos. Gracias a esta modificación se obtienen parámetros que varían en el tiempo y en el espacio lo que permite capturar la variabilidad espacial de la precipitación. Al igual que en el modelo de ocurrencia, la inclusión del total trimestral es necesaria si se desean simular series condicionadas.
Modelación: Clima local: Temperatura
Para el caso de la temperatura se utiliza una metodología similar a la utilizada para la precipitación, basada en Kleiber et al. (2013). A partir de los datos observados de temperatura se ajustan dos modelos: uno de máxima y otro de mínima. Ambos modelos utilizan las mismas variables para realizar el ajuste. El modelo puede ser expresado de la siguiente manera:
\[ X_{s,t} = (s, O_{s,t}, O_{s,t-1}, f(T_{x_{(s, t-1)}}, T_{n_{(s, t-1)}}), f(doy(t)), f(SX(t), SN(t)), f(lon, lat)) \] donde, \(X_{s,t}\) corresponde a la temperatura máxima o mínima en un sitio y momento determinado; \(s\) corresponde al efecto del sitio s (ordenada al origen); \(O_{s,t}\) corresponde a la ocurrencia de lluvia en sitio y día determinado; \(O_{s,t-1}\) corresponde a la ocurrencia del día previo; \(f(T_{x_{(s, t-1)}}, T_{n_{(s, t-1)}})\) corresponde a una función suavizada de la interacción entre la temperatura máxima y mínima del día previo; \(f(doy(t))\) corresponde a una función cíclica de los días del año para considerar el efecto de la estacionalidad sobre la temperatura; El término \(f(SX(t), SN(t))\) corresponde a la interacción entre las medias estacionales de temperatura máxima y mínima y, como se explicó en la sección anterior, se incluyen para crear modelos condicionados y \(f(lon, lat)\) corresponde a las coordenadas del punto s (solo se utiliza en el modelo espacial). La ocurrencia de lluvia es muy importante ya que en general los días lluviosos tienen temperaturas más bajas que los secos, sobre todo en verano, por lo que debe ser incluido en el ajuste. La Figura 6 es un ejemplo de esta influencia en Junín, en donde se muestra la amplitud térmica para cada mes del año en función del tipo de día, seco o lluvioso.
Modelación: Tiempo local: Para una estación
Precipitación
El tiempo local de la ocurrencia de precipitación se modela a través de los residuos de la regresión probit. Por definición estos residuos tienen una distribución \(X \sim \mathcal{N}(0,\,1)\,\). Por lo tanto, generar valores para el tiempo local es muy sencillo, solo se necesita de una función gaussiana que genere números aleatorios. La Figura muestra un ejemplo de la generación de clima local y tiempo local para una estación meteorológica de Argentina.
Figura 6. Tiempo local: precipitación
Esta Figura muestra sólo un año para mejorar la visualización pero su interpretación es válida para toda la longitud de la serie. La línea azul corresponde al clima local modelado a través del GAM y la naranja al ruido aleatorio creado a partir de una distribución \(X \sim \mathcal{N}(0,\,1)\,\). La sumatoria de ambos componentes determinarán si el día es lluvioso o no. Si la suma es positiva, lloverá, caso contrario será un día seco. Se puede observar en la línea azul un patrón estacional debido al régimen de precipitación tipo monzónico de esta región con picos de más días lluviosos durante el verano mientras que en invierno disminuyen marcadamente. Esta misma serie temporal de números aleatorios serán la base para la generación de los montos de precipitación.
Temperatura
El tiempo local de las temperaturas máxima y mínima se modela de una manera diferente. En este caso se toman los residuos de cada uno de los GAMs, es decir, la diferencia entre el valor ajustado por el modelo y el valor observado de temperatura. Además, como la temperatura está fuertemente influenciada por el tipo de día, días lluviosos tienden a tener una menor amplitud térmica que los días secos, los residuos se separan en función del tipo de día. Para capturar mejor el patrón estacional de la temperatura, los residuos se agrupan por mes y se ajusta un modelo bivariado que contemple los residuos de temperaturas máxima y mínima. El uso de un modelos bivariado es una alternativa muy interesante para mantener la consistencia entre ambas variables, es decir, que la temperatura máxima no supere a la mínima. La siguiente Figura es un ejemplo de las series de tiempo local para una localidad de Argentina.
Figura 7. Tiempo local: temperatura
La figura esta dividida en dos paneles, el superior muestra el tiempo local para los días lluviosos y el inferior para los días secos. La línea naranja corresponde a los valores de tiempo local para temperatura máxima y los azules para temperatura mínima. Si bien ambas series difieren en magnitud y variabilidad, las dos tienden a variar conjuntamente.
Modelación: Tiempo local: Para una grilla
Para generar datos sobre una grilla regular o en puntos donde no hay datos para ajustar los modelos se debe utilizar el modelo espacial. La generación del tiempo local en el espacio es conceptualmente idéntico a la generación sobre estaciones meteorológicas sólo que utiliza campos gaussianos para incluir la dimensión espcial. Los campos gaussianos se generan con el paquete RandomFields y capturan la variabilidad espacial de cada una de las variables a partir de su variograma.
Precipitación
Para la precipitación se utiliza un modelo exponencial que utiliza como parámetros el variograma ajustado a partir de los datos observados usando máxima verosimilitud. Este modelo permite simular campos con una muy buena consistencia espacial. La ocurrencia de precipitación no ocurre de manera aislada en una región sino que, en general, un evento lluvioso abarca una importante superficie. La siguiente Figura es una ejemplo de los campos aleatorios generados sobre una grilla regular para una región de Argentina.
Figura 8. Tiempo local: precipitación sobre una grilla
La interpretación es análoga a la mostrada para una estación puntual. Los valores para cada píxel se suman a la componente climática y así se determina si el día será lluvioso o no. Este tipo de campos se generan para todos los días de la simulación y cada campo diario es independiente del anterior.
Temperatura
Al igual que para una estación, los campos gaussianos que se generan son bivariados. Se utiliza el modelo Bivariado de Whittle Matern incluido en el paquete RandomFields. La siguiente Figura es un ejemplo para un día de una realización para grilla regular sobre Argentina.
Figura 8. Tiempo local: temperatura sobre una grilla
Ambos campos se generan simultáneamente para temperaturas máxima y minima y discriminando entre días secos y lluviosos. Esto quiere decir que en el proceso de creación de los campos se generan cuatro capas diferentes que luego se combinan en función del tipo de día de cada píxel resultando en dos campos integradores. Los valores para cada pi’xel se sumen a la componente climática y así se obtiene el valor final para cada día de la simulación.
Aplicación
En esta sección se mostrarán ejemplos de aplicación del generador para generar distintos tipos de series explicando las funciones necesarias y cada uno de los parámetros.
Instalar paquetes necesarios
El primer paso es comprobar que todos los paquetes necesarios estén instalados y si no es así, descargarlos e intalarlos.
## Loading required package: gamwgenCreación de directorios
El paquete tiene precargados algunos ejemplos de aplicación con datos reales de estaciones meteorológicas de la red del SISSA.
Para poder seguir este tutorial se deben crear directorios donde se guardarán los datos de entrada y salida.
- /input_data: aquí se guardarán los datos meteorológicos y los metadatos de las estaciones
- /output_data: aquí se guardarán los resultados de la simulación
Si estos directorios no existen, se crearán.
Generación de series sintéticas para una sola estación meteorológica
Este primer ejemplo consiste en la generación de series sintéticas para estaciones meteorológicas. Es decir, se generarán series para las mismas estaciones que fueron utilizadas en el ajuste de los distintos modelos. Por lo tanto, no se considera la componente espacial. El ejemplo está divido en dos, en una primera parte se ajustará el modelo para generar series estacionarias y luego, en una segunda parte, se incluirán coviariables estacionales para producir series pseudohistóricas. Para ambos ejemplos los datos son usados serán los mismos.
Crear archivos de entrada
El primer paso consiste en generar los set de datos de entrada que se descargaron al momento de instalar el paquete del generador estocástico. Estos datos son sólo a título demostrativo, si el usuario desea correr el modelo con sus propios datos deberá cambiar los objetos que se generarán en esta sección por los suyos y colocarlos en la carpeta input_data.
Los archivos necesarios son:
- stations.csv
- climate.csv
Los datos meteorológicos se dividen en dos archivos separados: stations.csv y climate.csv. Los nombres de los mismos no deben ser necesariamente iguales a los usados aquí.
Los metadatos de las estaciones se alojan en el archivo stations.csv. Este archivo contiene la información de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:
- station_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integer
- latitude: latitud en grados decimales. La variable debe ser de tipo double
- longitude: longitud en grados decimales. La variable debe ser de tipo double
La tabla puede tener más variables pero sólo se necesitan las anteriores.
A continuación se muestran la primera fila del dataset y los tipo de datos de cada una de las variables.
| x | y | station_id | nombre | lat_dec | lon_dec | elev | pais_id |
|---|---|---|---|---|---|---|---|
| 5001614 | 6256841 | 87448 | Villa Reynolds Aero | -33.7181 | -65.3737 | 486 | AR |
El objeto stations debe ser convertido de tibble a sf. El sistema de referencia espacial debe ser planar. No es necesario un sistema de referencia espacial en particular, solamente las coordenadas deben estar expresadas en metros.
# Convertimos el objeto stations a sf y se convierte su proyección de WGS 1984 a
# POSGAR Argentina Faja 5.
stations %<>%
sf::st_as_sf(coords = c("lon_dec", "lat_dec"), crs = 4326) %>%
sf::st_transform(crs = 22185)La información climática se aloja en el archivo climate.csv. Este archivo contiene los datos de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:
date: fecha del dato. La variable debe ser de tipo datestation_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integerprcp: datos diarios de precipitación La variable debe ser de tipo doubletmax: datos diarios de temperatura máxima. La variable debe ser de tipo doubletmin: datos diarios de temperatura mínima. La variable debe ser de tipo double
A continuación se muestran las primeras cinco filas del dataset y los tipos de datos de cada una de las variables.
| date | station_id | tmax | tmin | prcp |
|---|---|---|---|---|
| 1961-01-01 | 87448 | 37.4 | 13.5 | 0.6 |
| 1961-01-02 | 87448 | 27.4 | 14.3 | 23.9 |
| 1961-01-03 | 87448 | 26.6 | 13.5 | 0.0 |
| 1961-01-04 | 87448 | 31.0 | 11.7 | 6.0 |
| 1961-01-05 | 87448 | 27.0 | 14.1 | 0.0 |
| 1961-01-06 | 87448 | 26.3 | 11.3 | 0.0 |
| 1961-01-07 | 87448 | 34.1 | 12.0 | 6.7 |
| 1961-01-08 | 87448 | 32.8 | 15.9 | 0.0 |
| 1961-01-09 | 87448 | 37.6 | 16.1 | 0.0 |
| 1961-01-10 | 87448 | 26.9 | 4.6 | 0.0 |
Los nombres de las variables son importantes y deben ser siempre los mismos ya que el modelo las reconocerá a partir de los mismos. Los nombres deben ser los siguientes:
date: corresponde a la fecha del día en formatoDate. El formato de la fecha para facilitar el reconocimiento por parte de R es “YYYY-MM-DD”, es decir, el año expresado con cuatro dígitos y luego dos dígitos para el mes y dos para el día.station_id: Identificador unívoco de cada una de las estaciones. Debe ser un número entero.tmax: temperatura máxima diaria expresada en °C.tmin: temperatura mínima diaria expresada en °C.prcp: precipitación diaria expresada en mm.
El orden de las variables no es importante pero, como se mencionó, si se deben respetar los nombres de cada una. En el caso de faltantes, no se utiliza ningún valor específico para los NAs, sólo se debe dejar ese valor vacío. Este archivo tiene un formato largo, es decir, las estaciones se deben colocar una debajo de la otra.
La generación de series sobre estaciones requiere de dos funciones básicas local_fit y local_simulate. Independientemente de si se incluyen totales trimestrales en el modelo, siempre se utilizan esas dos funciones.
Series sintéticas estacionarias
Ajuste de los modelos
A continuación se mostrará como ajustar los cuatro modelos estadísticos para una sola estación meteorológica para generar series estacionarias donde cada realización es complemente independiente.
El ajuste del modelo local (en un punto) necesita de dos funciones: en una se define la configuración general del modelo y con la segunda se corre el modelo propiamente dicho.
Primero se crea un objeto con el control para el ajuste del simulador. Los argumentos son:
prcp_occurrence_threshold: umbral de precipitación para un día lluvioso. La OMM recomienda un umbral de 0.1 mm para considerar un día como lluvioso.avbl_cores: cantidad de núcleos disponibles para la paralelización.planar_crs_in_metric_coords: sistema de coordenadas planar.
control_fit <-gamwgen::local_fit_control(
prcp_occurrence_threshold = 0.1, # Umbral para la definición de días húmedos
avbl_cores = 6, # Cantidad de núcleos disponibles
planar_crs_in_metric_coords = 22185) # Sistema de referencia espacial (en metros)Luego se corre el ajuste para la estación meteorológica con la función local_calibrate. Los argumentos de la función son:
climate: datos meteorológicos observados para la estaciónstations: metadatos de las estaciones meteorológicasseasonal_covariates: datos agregados trimestrales. Si es NULL el ajuste será sin covariables y las series generadas serán estacionarias.control: objeto de controlverbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza el ajuste de los cuatro modelos para cada una de
# las estaciones. En este caso, por cuestiones de tiempo a cargar un objeto ya precalculado.
# Si el usuario desea correrlo deberá ver la nota anterior.
gamgen_fit <- gamwgen::local_calibrate(climate = climate, # Registro histórico de variables meteorológicas
stations = stations, # Estaciones meteorológicas
seasonal_covariates = NULL, # Totales trimestrales de precipitación
control = control_fit, # Objeto de control
verbose = FALSE) # Impresión de mensajes en la consola.Para esta demostración, cargamos el objeto con el ajuste del modelo ya realizado.
# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('input_data/local/fit_local_unconditional.RData')) {
fs::file_copy(system.file('/autorun/local', "fit_local_unconditional.RData", package = "gamwgen"),
new_path = 'input_data/local/fit_local_unconditional.RData')
}
# Cargamos el archivo recientemente creado
load('input_data/local/fit_local_unconditional.RData')
# Clase del objeto con el ajuste del generador
class(gamgen_fit)## [1] "gamwgen"## [1] "control" "stations" "climate"
## [4] "crs_used_to_fit" "start_climatology" "fitted_models"
## [7] "models_data" "models_residuals" "statistics_threshold"
## [10] "exec_times"Dentro del objeto se guardan todo lo necesario para la simulación así como información accesoria.
control: copia de la configuración usada para calibrar el generadorstations: estaciones meteorológicas utilizadas para la calibraciónclimate: datos climáticos de cada uno de las estacionesseasonal_covariates: series temporales de totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima.crs_used_to_fit: sistema de referencia espacial usado para proyectarstart_climatology: climatología diaria de cada una de las variables de entrada.fitted_models: modelos ajustados, uno para cada variable: temperaturas máxima y mínima y ocurrencia y montos de precipitación.models_data: datos usados efectivamente usados para ajustar los modelos (sin NAs)models_residuals: residuos de cada uno de los modelos. Es decir, la diferencia entre el valor ajustado por el modelo (clima local) y el valor observado en el día istatistics_threshold: umbrales de amplitud térmica diaria por mes. Si la amplitud simulada está fuera de este rango, se repetirá la simulación para ese día a los fines de mantener la consistencia entre variablesexec_times: tiempo de ejecución de cada una de las etapas del ajuste
Cada uno de los GAMs ajustados se almacenan en el objeto gamgen_fit y pueden ser evaluados con la función summary().
##
## Family: gaussian
## Link function: identity
##
## Formula:
## tmax ~ s(tmax_prev, tmin_prev, k = 50) + s(prcp_occ, bs = "re") +
## s(prcp_occ_prev, bs = "re") + s(doy, bs = "cc", k = 30)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.61343 0.03214 796.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(tmax_prev,tmin_prev) 30.2469 38.71 247.1 <2e-16 ***
## s(prcp_occ) 0.9989 1.00 1005.8 <2e-16 ***
## s(prcp_occ_prev) 0.9995 1.00 2300.2 <2e-16 ***
## s(doy) 12.8918 28.00 158.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.71 Deviance explained = 71.1%
## fREML = 58955 Scale est. = 14.138 n = 21466La función summary permite analizar los resultados del ajuste de cada uno de los modelos. Para el caso del modelo de temperatura máxima podemos ver la fórmula del GAM en la parte superior bajo el apartado Formula y la significancia de cada uno de los términos del modelo en la tabla inmediatamente inferior. Se puede observar que todos los términos son altamente significativos. También se incluyen como pruebas de bondad del ajuste el porcentaje de la varianza explicada por el modelo y el valor de R-ajustado.
El paquete gratia (Simpson, Gavin (2018)) es muy útil para la visualización de los ajustes de manera gráfica. Esta función toma los diagnósticos por defecto de la función gam.check del paquete mgcv (Wood, Simon (2001)) y los convierte en gráficos de la librería ggplot2 (Wickham, Hadley(2011)) que son visualmente muy atractivos. La figura está dividida en cuatro paneles, cada uno mostrando un diagnóstico diferente. En el panel superior izquierdo se muestra un Q-Q Plot de los residuos. Si el modelo ha ajustado satisfactoriamente los datos, los puntos deberían estar sobre la recta 1:1 demarcada en rojo. El panel superior derecho muestra los residuos vs el predictor lineal. El objetivo de este diagnóstico es que los residuos se distribuyan al azar y que no tengan un patrón claro. La presencia de patrones en los residuos indicaría que el modelo no ha explicado algun componente importante de la variabilidad de los datos. En el panel inferior izquierdo se muestra un histograma de los residuos siendo deseable que lo mismos tengan una distribución proxima a la Normal. El último gráfico en el panel inferior derecho muestra una gráfico de dispersión entre los valores ajustados y observados.
Estos diagnóstico exploratorios puede aplicarse a cada uno de los modelos ajustados por la función local_calibrate.
Con la creación del objeto gamgen_fit finaliza el proceso de ajuste y ya se pueden generar series sintéticas para la o las estaciones usadas para calibrar el generador.
Generación de series
La generación de series sigue la misma estructura anterior, una función para configurar la generación y otra que realiza la generación propiamente dicha.
Los argumentos de la función de control son:
nsim: cantidad de simulaciones a realizar. Se debe ingresar un valor entero mayor o igual a 1.seed: semilla. Se debe ingresar cualquier numero entero. No es necesario recordarlo porque se guarda junto a los resultados.avbl_cores: cantidad de núcleos disponibles para la paralelización.use_spatially_correlated_noise: utilizar la generación estocástica espacialmente correlacionada. Esta opción sólo es válida si en el ajuste y en la simulación se usaron más de cinco estaciones meteorológicas diferentes. Con un menor número no es posible calcular los variogramas necesarios para la generación de los campos aleatorios. Se debe introducir un boolean (TRUE or FALSE).use_temporary_files_to_save_ram: si se simulan muchas realizaciones o los recursos informáticos son escasos, esta opción permite guardar los resultados de cada una de las realizaciones en el disco liberando memoria RAM que quedará disponible para generar nuevas simulaciones. Al finalizar la generación todos los archivos se combinan en uno único. Se debe introducir un boolean (TRUE or FALSE).use_temporary_files_to_save_ram: esta opción permite eliminar los archivos temporales creados para ahorrar RAM luego de terminar la generación de todas las simulaciones. Se debe introducir un boolean (TRUE or FALSE).
control_sim <- gamwgen::local_simulation_control(
nsim = 10, # Cantidad de simulaciones a realizar
seed = 1234, # Semilla para que los resultados sean reproducibles
avbl_cores = 6, # Cantidad de núcleos disponibles a utilizar
use_spatially_correlated_noise = FALSE, # Usar modelo de ruido espacialmente correlacionado
use_temporary_files_to_save_ram = FALSE, # Guardar resultados intermedios para ahorrar RAM
remove_temp_files_used_to_save_ram = TRUE) # Borrar los resultados intermedios creados anteriormenteLuego se procede a la simulación de datos meteorológicos. Los argumentos de la función de simulación son:
model: objeto con el resultado de la funciónlocal_calibrate()simulation_locations: objeto tiposfcon la ubicación de las estaciones a simular. Las estaciones usadas deben haber sido incluidas en el proceso de ajuste. No es necesario que todas estén presentes, se pueden generar series solo sobre algunas de ellas.start_date: fecha de comienzo de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de comienzo es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, ni tampoco posterior. Se debe introducir una fecha en formato dateend_date: fecha de fin de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de fin es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, tampoco puede ser posterior. Se debe introducir una fecha en formato datecontrol: objeto de control creado con la funcióncontrol_sim().output_folder: ruta al directorio donde se guardarán los resultados, tanto finales como intermedios.output_filename: nombre del archivo de salida. Para facilitar la interoperabilidad, el archivo generado es un archivo de texto en formato separado por comas (.csv)seasonal_covariates: datos agregados trimestrales. Si el ajuste se realizó con covariables, la generación también debe realizarse con ellas. Caso contrario se producirá un error. Se debe introducir un data frame con los valores agregados para las tres variables (precipitación y temperaturas máxima y mínima) pero no necesariamente deben ser los mismos a los utilizados en el ajuste. Si se desean simular tendencias de algún tipo, ya sea de un modelo de cambio climático o arbitrarias, se deben perturbar estas variables trimestrales e introducirlas aquí. Estas series si deben tener la misma longitud que el período a generarverbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza la generación de series para cada una de las estaciones.
# En este caso, por cuestiones de tiempo, vamos a cargar un objeto con los resultados de la simulación
simulated_climate <- gamwgen::local_simulation(model = gamgen_fit, # Objeto con los resultados del ajuste
simulation_locations = stations, # Estaciones para las cuales simular
start_date = as.Date('2019-01-01'), # Fecha de comienzo de las simulaciones
end_date = as.Date('2019-12-31'), # Fecha de fin de las simulaciones
control = control_sim, # Objeto con la configuración
output_folder = getwd(), # Directorio donde se guardarán los resultados
output_filename = 'simulations.csv', # Nombre del archivo de salida
seasonal_covariates = NULL, # Covariables estacionales
verbose = FALSE) # Impresión de mensajes en la consolaEsta función produce dos tipos de resultados: una lista que permanece en el ambiente de R y los datos generados que son guardados como .csv en el directorio indicado precedentemente.
# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('output_data/simulated_local_unconditional.RData')) {
fs::file_copy(system.file('/autorun/local', "simulated_local_unconditional.RData", package = "gamwgen"),
new_path = 'output_data/simulated_local_unconditional.RData')
}
# Cargamos el archivo recientemente creado
load('output_data/simulated_local_unconditional.RData')
# Clase del objeto con el ajuste del generador
class(simulated_climate)## [1] "list" "gamwgen.climate"## [1] "nsim"
## [2] "seed"
## [3] "realizations_seeds"
## [4] "simulation_points"
## [5] "output_file_with_results"
## [6] "output_file_fomart"
## [7] "rdata_file_with_fitted_stations_and_climate"
## [8] "exec_times"La lista contiene los siguientes objetos:
nsim: cantidad de realizaciones.seed: semilla general para toda la generación. Corresponde a la que se incluye en la función de control.realization_seeds: semillas para cada una de las realizaciones. Esto permite replicar los resultados.simulation_points: puntos donde se generaron las series sintéticas.output_fil_with_results: nombre del archivo con los resultados.output_file_format: tipo de archivo de salida, en este caso .csv.rdata_file_with_fitted_stations_and_climate: archivo.RDatacon los datos meteorológicos observados que fueron utilizados en el ajuste. También se incluyen los metadatos de cada uno de esos puntos.exec_times: tiempo de ejecución del la generación.
Ahora veremos el formato del archivo de salida que contiene las series sintéticas.
Para desmenuzar la generación del tiempo local se pueden correr algunas de las funciones que forman parte de local_simulation. Por ejemplo, del objeto gamgen_fit se extraen los residuos y con la función gamwgen:::generate_residuals_statistics se calculan los parámetros.
En la tabla anterior se muestran los parámetros necesarios para generar el tiempo local de las temperaturas máxima y mínima.
gen_noise_params <- gamwgen:::generate_residuals_statistics(
models_residuals = gamgen_fit$models_residuals)
#knitr::kable(gen_noise_params)
rmarkdown::paged_table(gen_noise_params)station_id: número unívoco que identifica a cada estación meteorológica.type: tipo de día lluvioso (Wet) o seco (Dry).month: número de mes para los que se calculan los parámetrossd.tmax_residuals: desvío estándar de los residuos del modelo de temperatura máxima.sd.tmin_residuals: desvío estándar de los residuos del modelo de temperatura mínima.mean.tmax_residuals: media de los residuos del modelo de temperatura máxima.mean.tmin_residuals: media de los residuos del modelo de temperatura mínima *cov.residuals: covarianza de los residuos.var.tmax_residuals: covarianza de los residuos del modelo de temperatura máxima.var.tmin_residuals: covarianza de los residuos del modelo de temperatura mínima.
Losa parámetros para cada uno de los meses permiten generar valores de tiempo local para las dos temperaturas a partir de una distribución normal multivariada.
A continuación se muestra un ejemplo para el año 2019 sólo para los días lluviosos.
fechas <- data.frame(date = seq(as.Date('2019-01-01'), as.Date('2019-12-31'), 'days')) %>%
dplyr::mutate(dia = lubridate::day(date),
mes = lubridate::month(date))
# Ejemplo de generación de ruido para el mes de enero
tmax_dry <- purrr::map2_dfr(
.x = fechas$dia,
.y = fechas$mes,
.f = function(dia, mes) {
result_tmax_dry <- control_sim$temperature_noise_generating_function(simulation_points = stations %>%
dplyr::filter(., station_id == '87448'),
gen_noise_params = gen_noise_params,
month_number = mes,
selector = 'tmax_dry',
seed = NULL)
result_tmax_dry <- result_tmax_dry %>%
dplyr::mutate(dia = dia,
mes = mes)
}
) %>%
sf::st_set_geometry(NULL) %>%
dplyr::mutate(date = as.Date(paste0('2019-', mes, '-', dia))) %>%
dplyr::select(-mes, -dia)
ggplot2::ggplot(data = tmax_dry %>%
tidyr::gather(residuo, valor, -date),
ggplot2::aes(x = date, y = valor, color = residuo)) +
ggplot2::scale_y_continuous(limits = c(-15, 15),
breaks = seq(-15, 15, 3),
name = 'Tiempo local') +
ggplot2::scale_x_date(name = 'Días') +
ggplot2::scale_color_manual(values=c("DarkOrange", "Steelblue"),
labels = c("Temp. Máx.", "Temp. Min.")) +
ggplot2::geom_line() +
ggplot2::theme_bw() +
ggplot2::theme(legend.position="bottom",
legend.title = ggplot2::element_blank())
La líonea naranja corresponde a la serie para temperaturas máximas y la azul para temperaturas mínimas.
# Se carga el set de datos simulados
simulated_climate <- readr::read_csv(here::here('output_data/local/simulated_local_unconditional.csv')) ##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## realization = col_double(),
## station_id = col_double(),
## point_id = col_double(),
## longitude = col_double(),
## latitude = col_double(),
## date = col_date(format = ""),
## tmax = col_double(),
## tmin = col_double(),
## prcp = col_double()
## )| realization | station_id | point_id | longitude | latitude | date | tmax | tmin | prcp |
|---|---|---|---|---|---|---|---|---|
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-01 | 27.19266 | 5.464899 | 0 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-02 | 30.57296 | 10.266542 | 0 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-03 | 33.36446 | 15.180860 | 0 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-04 | 34.05627 | 15.473379 | 0 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-05 | 32.69240 | 15.270145 | 0 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-06 | 30.56683 | 14.797530 | 0 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-07 | 29.32913 | 12.624675 | 0 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-08 | 34.09988 | 13.816714 | 0 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-09 | 35.03660 | 19.596069 | 0 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2019-01-10 | 37.52834 | 19.323512 | 0 |
El resultado de la generación es un archivo .csv que contiene la siguiente información:
realization: número de realización. Es un valor entero entre 1 y la cantidad de realizaciones definida por el usuario.station_id: número unívoco de identificación de la estación meteorológica o del punto arbitrario.date: fechas de cada uno de los días de la simulación.tmax: valores de temperatura máxima generada expresada en °C.tmin: valores de temperatura mínima generada expresada en °C.prcp: valores de precipitación diaria generada expresada en mm.
La siguiente Figura muestra un ejemplo de las series de temperaturas máximas y mínimas generadas.
# Grafico de temperatura máxima
tmax_plot <- ggplot2::ggplot() +
ggplot2::geom_line(data = simulated_climate,
ggplot2::aes(x = date, y = tmax, color = realization),
alpha = 0.5, color = 'DarkOrange') +
ggplot2::geom_line(data = climate %>%
dplyr::filter(date > as.Date('2019-01-01')),
ggplot2::aes(x = date, y = tmax)) +
ggplot2::labs(x = 'Fecha', y = 'Temperatura máxima [°C]') +
ggplot2::theme_bw() +
ggplot2::theme(legend.position = "none")
# Grafico de temperatura mínima
tmin_plot <- ggplot2::ggplot() +
ggplot2::geom_line(data = simulated_climate,
ggplot2::aes(x = date, y = tmin, color = realization),
alpha = 0.5, color = 'Steelblue') +
ggplot2::geom_line(data = climate %>%
dplyr::filter(date > as.Date('2019-01-01')),
ggplot2::aes(x = date, y = tmin)) +
ggplot2::labs(x = 'Fecha', y = 'Temperatura mínima [°C]') +
ggplot2::theme_bw() +
ggplot2::theme(legend.position = "none")
cowplot::ggdraw() +
cowplot::draw_plot(tmax_plot, x = 0, y = 0.5, width = 1, height = .5) +
cowplot::draw_plot(tmin_plot, x = 0, y = 0, width = 1, height = .5)
En el panel superior se muestra la temperatura máxima observada (negra) y las temperaturas sintéticas (naranja) para el año 2019. En el infeior se muestra la temperatura mínima diaria observada (negra) y cada una de las realizaciones (azul).
Series sintéticas pseudohistóricas
Ajuste de los modelos
A continuación se mostrará como ajustar los cuatro modelos estadísticos para una sola estación meteorológica para generar series pseudohistóricas donde cada realización copia las variaciones de baja frecuencia e la serie observada. El ajuste del modelo local (en un punto) necesita de dos funciones: en una se define la configuración general del modelo y con la segunda se corre el modelo propiamente dicho.
Primero se crea un objeto con el control para el ajuste del simulador. Los argumentos son:
prcp_occurrence_threshold: umbral de precipitación para un día lluvioso. La OMM recomienda un umbral de 0.1 mm para considerar un día como lluvioso.avbl_cores: cantidad de núcleos disponibles para la paralelización.planar_crs_in_metric_coords: sistema de coordenadas planar.
control_fit <-gamwgen::local_fit_control(
prcp_occurrence_threshold = 0.1, # Umbral para la definición de días húmedos
avbl_cores = 6, # Cantidad de núcleos disponibles
planar_crs_in_metric_coords = 22185) # Sistema de referencia espacial (en metros)Al tratarse de un modelo que ajusta condicionado por la variabilidad de baja frecuencia preexistente en los datos observados es necesario la agregación de las variables diarias en totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima. Esta operación puede realizarse con la función summarise_seasonal_climate incluida en el paquete. Esta función, además de agregar los datos, permite la imputación de faltantes. Se toleran una cierta cantidad que puede ser determinada por el usuario. El método de imputación utilizado es el imputePCA() de la librería missMDA .
# Agregación de valores diarios
seasonal_covariates <- gamwgen::summarise_seasonal_climate(climate, umbral_faltantes = 0.2)
# Se muestran las primeras cinco filas
knitr::kable(seasonal_covariates[1:10,])| station_id | year | season | seasonal_prcp | seasonal_tmax | seasonal_tmin |
|---|---|---|---|---|---|
| 87448 | 1961 | 1 | 273.7 | 30.98764 | 14.5134831 |
| 87448 | 1961 | 2 | 236.0 | 24.47473 | 8.4670330 |
| 87448 | 1961 | 3 | 11.3 | 19.04565 | 1.3456522 |
| 87448 | 1961 | 4 | 234.5 | 24.66235 | 7.7611765 |
| 87448 | 1962 | 1 | 402.4 | 29.92333 | 14.2455556 |
| 87448 | 1962 | 2 | 187.2 | 24.30440 | 7.9868132 |
| 87448 | 1962 | 3 | 50.9 | 17.21957 | -0.8978261 |
| 87448 | 1962 | 4 | 179.7 | 25.39560 | 7.8472527 |
| 87448 | 1963 | 1 | 306.8 | 29.09101 | 14.1404494 |
| 87448 | 1963 | 2 | 82.0 | 25.46196 | 8.5695652 |
Cabe mencionar que con esta función las valores se agregan por trimestre considerando la siguiente definición:
- Verano: Diciembre, Enero y Febrero
- Otoño: Marzo, Abril y Mayo
- Invierno: Junio, Julio y Agosto
- Primavera: Septiembre, Octubre y Noviembre
Los valores también se podrían agregar siguiendo otra definición de estaciones pero en ese caso, el usuario debería hacerlo por su cuenta. Algunas funciones útiles para hacerlo son las disponibles en el paquete lubridate como quarter() que permite definir el mes de comienzo de los trimestres. Para estas variables los nombres también son importantes por lo que deben respetarse los mostrados anteriormente.
La siguiente Figura muestra los totales trimestrales de precipitación para la estación meteorológica del ejemplo.
ggplot2::ggplot(data = seasonal_covariates %>%
dplyr::mutate(season = factor(season,
labels = c('Verano', 'Otoño', 'Invierno', 'Primavera'))),
ggplot2::aes(x = year, y = seasonal_prcp, group = 1)) +
ggplot2::geom_line() +
ggplot2::geom_smooth(se = FALSE, span = 0.5) +
ggplot2::facet_wrap(~season, scales = 'free') +
ggplot2::labs(x = 'Años', y = 'Precipitación acumulada [mm]') +
ggplot2::theme_bw()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Cada uno de los paneles muestra la precipitación acumulada por trimestre junto a una regresión local (loess) en azul. Al incluir estos totales en el modelo, las series generadas tenderán a seguir la variabilidad observada, es decir, años secos generarán realización con valores nferiores a la media y viceversa.
Luego se corre el ajuste para la estación meteorológica con la función local_calibrate. Los argumentos de la función son:
climate: datos meteorológicos observados para la estaciónstations: metadatos de las estaciones meteorológicasseasonal_covariates: datos agregados trimestrales. Si es NULL el ajuste será sin covariables y las series generadas serán estacionarias.control: objeto de controlverbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza el ajuste de los cuatro modelos para cada una de
# las estaciones. En este caso, por cuestiones de tiempo a cargar un objeto ya precalculado.
# Si el usuario desea correrlo deberá ver la nota anterior.
gamgen_fit <- gamwgen::local_calibrate(climate = climate, # Registro histórico de variables meteorológicas
stations = stations, # Estaciones meteorológicas
seasonal_covariates = seasonal_covariates, # Totales trimestrales de precipitación
control = control_fit, # Objeto de control
verbose = FALSE) # Impresión de mensajes en la consola.Para esta demostración, cargamos el objeto con el ajuste del modelo ya realizado.
# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('input_data/local/fit_local_conditional.RData')) {
fs::file_copy(system.file('/autorun/local', "fit_local_conditional.RData", package = "gamwgen"),
new_path = 'input_data/local/fit_local_conditional.RData')
}
# Cargamos el archivo recientemente creado
load('input_data/local/fit_local_conditional.RData')
# Clase del objeto con el ajuste del generador
class(gamgen_fit)## [1] "gamwgen"## [1] "control" "stations" "climate"
## [4] "seasonal_covariates" "crs_used_to_fit" "start_climatology"
## [7] "fitted_models" "models_data" "models_residuals"
## [10] "statistics_threshold" "exec_times"Dentro del objeto se guardan todo lo necesario para la simulación así como información accesoria.
control: copia de la configuración usada para calibrar el generadorstations: estaciones meteorológicas utilizadas para la calibraciónclimate: datos climáticos de cada uno de las estacionesseasonal_covariates: series temporales de totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima.crs_used_to_fit: sistema de referencia espacial usado para proyectarstart_climatology: climatología diaria de cada una de las variables de entrada.fitted_models: modelos ajustados, uno para cada variable: temperaturas máxima y mínima y ocurrencia y montos de precipitación.models_data: datos usados efectivamente usados para ajustar los modelos (sin NAs)models_residuals: residuos de cada uno de los modelos. Es decir, la diferencia entre el valor ajustado por el modelo (clima local) y el valor observado en el día istatistics_threshold: umbrales de amplitud térmica diaria por mes. Si la amplitud simulada está fuera de este rango, se repetirá la simulación para ese día a los fines de mantener la consistencia entre variablesexec_times: tiempo de ejecución de cada una de las etapas del ajuste
Cada uno de los GAMs ajustados se almacenan en el objeto gamgen_fit y pueden ser evaluados con la función summary().
##
## Family: gaussian
## Link function: identity
##
## Formula:
## tmax ~ s(tmax_prev, tmin_prev, k = 50) + s(prcp_occ, bs = "re") +
## s(prcp_occ_prev, bs = "re") + s(doy, bs = "cc", k = 30) +
## s(SX1, SN1, k = 20) + s(SX2, SN2, k = 20) + s(SX3, SN3, k = 20) +
## s(SX4, SN4, k = 20)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.62116 0.03189 803.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(tmax_prev,tmin_prev) 29.9069 38.405 225.74 <2e-16 ***
## s(prcp_occ) 0.9989 1.000 930.30 <2e-16 ***
## s(prcp_occ_prev) 0.9995 1.000 2248.11 <2e-16 ***
## s(doy) 12.6336 28.000 73.12 <2e-16 ***
## s(SX1,SN1) 3.0997 3.689 55.32 <2e-16 ***
## s(SX2,SN2) 2.0001 2.000 89.25 <2e-16 ***
## s(SX3,SN3) 4.4053 5.853 35.79 <2e-16 ***
## s(SX4,SN4) 2.0001 2.000 86.85 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.713 Deviance explained = 71.4%
## fREML = 58829 Scale est. = 13.964 n = 21466La función summary permite analizar los resultados del ajuste de cada uno de los modelos. Para el caso del modelo de temperatura máxima podemos ver la fórmula del GAM en la parte superior bajo el apartado Formula y la significancia de cada uno de los términos del modelo en la tabla inmediatamente inferior. En este configuracion, al incluirle variables estacionales, se incoporan nuevos términos el modelo como SX1, SX2, SX3, SX4 y SN1, SN2, SN3, SN4, que corresponden a variables dummy de las medias trimestrales de temperatura máxima y mínima, respectivamente. Se puede observar que todos los términos son altamente significativos. También se incluyen como pruebas de bondad del ajuste el porcentaje de la varianza explicada por el modelo y el valor de R-ajustado.
El paquete gratia (Simpson, Gavin (2018)) es muy útil para la visualización de los ajustes de manera gráfica. Esta función toma los diagnósticos por defecto de la función gam.check del paquete mgcv (Wood, Simon (2001)) y los convierte en gráficos de la librería ggplot2 (Wickham, Hadley(2011)) que son visualmente muy atractivos. La figura está dividida en cuatro paneles, cada uno mostrando un diagnóstico diferente. En el panel superior izquierdo se muestra un Q-Q Plot de los residuos. Si el modelo ha ajustado satisfactoriamente los datos, los puntos deberían estar sobre la recta 1:1 demarcada en rojo. El panel superior derecho muestra los residuos vs el predictor lineal. El objetivo de este diagnóstico es que los residuos se distribuyan al azar y que no tengan un patrón claro. La presencia de patrones en los residuos indicaría que el modelo no ha explicado algun componente importante de la variabilidad de los datos. En el panel inferior izquierdo se muestra un histograma de los residuos siendo deseable que lo mismos tengan una distribución proxima a la Normal. El último gráfico en el panel inferior derecho muestra una gráfico de dispersión entre los valores ajustados y observados.
Estos diagnóstico exploratorios puede aplicarse a cada uno de los modelos ajustados por la función local_calibrate.
Con la creación del objeto gamgen_fit finaliza el proceso de ajuste y ya se pueden generar series sintéticas para la o las estaciones usadas para calibrar el generador.
Generación de series
La generación de series sigue la misma estructura anterior, una función para configurar la generación y otra que realiza la generación propiamente dicha.
Los argumentos de la función de control son:
nsim: cantidad de simulaciones a realizar. Se debe ingresar un valor entero mayor o igual a 1.seed: semilla. Se debe ingresar cualquier numero entero. No es necesario recordarlo porque se guarda junto a los resultados.avbl_cores: cantidad de núcleos disponibles para la paralelización.use_spatially_correlated_noise: utilizar la generación estocástica espacialmente correlacionada. Esta opción sólo es válida si en el ajuste y en la simulación se usaron más de cinco estaciones meteorológicas diferentes. Con un menor número no es posible calcular los variogramas necesarios para la generación de los campos aleatorios. Se debe introducir un boolean (TRUE or FALSE).use_temporary_files_to_save_ram: si se simulan muchas realizaciones o los recursos informáticos son escasos, esta opción permite guardar los resultados de cada una de las realizaciones en el disco liberando memoria RAM que quedará disponible para generar nuevas simulaciones. Al finalizar la generación todos los archivos se combinan en uno único. Se debe introducir un boolean (TRUE or FALSE).use_temporary_files_to_save_ram: esta opción permite eliminar los archivos temporales creados para ahorrar RAM luego de terminar la generación de todas las simulaciones. Se debe introducir un boolean (TRUE or FALSE).
control_sim <- gamwgen::local_simulation_control(
nsim = 10, # Cantidad de simulaciones a realizar
seed = 1234, # Semilla para que los resultados sean reproducibles
avbl_cores = 6, # Cantidad de núcleos disponibles a utilizar
use_spatially_correlated_noise = FALSE, # Usar modelo de ruido espacialmente correlacionado
use_temporary_files_to_save_ram = FALSE, # Guardar resultados intermedios para ahorrar RAM
remove_temp_files_used_to_save_ram = TRUE) # Borrar los resultados intermedios creados anteriormenteLuego se procede a la simulación de datos meteorológicos. Los argumentos de la función de simulación son:
model: objeto con el resultado de la funciónlocal_calibrate()simulation_locations: objeto tiposfcon la ubicación de las estaciones a simular. Las estaciones usadas deben haber sido incluidas en el proceso de ajuste. No es necesario que todas estén presentes, se pueden generar series solo sobre algunas de ellas.start_date: fecha de comienzo de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de comienzo es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, ni tampoco posterior. Se debe introducir una fecha en formato dateend_date: fecha de fin de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de fin es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, tampoco puede ser posterior. Se debe introducir una fecha en formato datecontrol: objeto de control creado con la funcióncontrol_sim().output_folder: ruta al directorio donde se guardarán los resultados, tanto finales como intermedios.output_filename: nombre del archivo de salida. Para facilitar la interoperabilidad, el archivo generado es un archivo de texto en formato separado por comas (.csv)seasonal_covariates: datos agregados trimestrales. Si el ajuste se realizó con covariables, la generación también debe realizarse con ellas. Caso contrario se producirá un error. Se debe introducir un data frame con los valores agregados para las tres variables (precipitación y temperaturas máxima y mínima) pero no necesariamente deben ser los mismos a los utilizados en el ajuste. Si se desean simular tendencias de algún tipo, ya sea de un modelo de cambio climático o arbitrarias, se deben perturbar estas variables trimestrales e introducirlas aquí. Estas series si deben tener la misma longitud que el período a generarverbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza la generación de series para cada una de las estaciones.
# En este caso, por cuestiones de tiempo, vamos a cargar un objeto con los resultados de la simulación
simulated_climate <- gamwgen::local_simulation(model = gamgen_fit, # Objeto con los resultados del ajuste
simulation_locations = stations, # Estaciones para las cuales simular
start_date = as.Date('2010-01-01'), # Fecha de comienzo de las simulaciones
end_date = as.Date('2019-12-31'), # Fecha de fin de las simulaciones
control = control_sim, # Objeto con la configuración
output_folder = getwd(), # Directorio donde se guardarán los resultados
output_filename = 'simulations.csv', # Nombre del archivo de salida
seasonal_covariates = seasonal_covariates, # Covariables estacionales
verbose = FALSE) # Impresión de mensajes en la consolaEsta función produce dos tipos de resultados: una lista que permanece en el ambiente de R y los datos generados que son guardados como .csv en el directorio indicado precedentemente.
# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('output_data/local/simulated_local_conditional.RData')) {
fs::file_copy(system.file('/autorun/local', "simulated_local_conditional.RData", package = "gamwgen"),
new_path = 'output_data/local/simulated_local_conditional.RData')
}
# Cargamos el archivo recientemente creado
load('output_data/local/simulated_local_conditional.RData')
# Clase del objeto con el ajuste del generador
class(simulated_climate)## [1] "list" "gamwgen.climate"## [1] "nsim"
## [2] "seed"
## [3] "realizations_seeds"
## [4] "simulation_points"
## [5] "output_file_with_results"
## [6] "output_file_fomart"
## [7] "rdata_file_with_fitted_stations_and_climate"
## [8] "exec_times"La lista contiene los siguientes objetos:
nsim: cantidad de realizaciones.seed: semilla general para toda la generación. Corresponde a la que se incluye en la función de control.realization_seeds: semillas para cada una de las realizaciones. Esto permite replicar los resultados.simulation_points: puntos donde se generaron las series sintéticas.output_fil_with_results: nombre del archivo con los resultados.output_file_format: tipo de archivo de salida, en este caso .csv.rdata_file_with_fitted_stations_and_climate: archivo.RDatacon los datos meteorológicos observados que fueron utilizados en el ajuste. También se incluyen los metadatos de cada uno de esos puntos.exec_times: tiempo de ejecución del la generación.
Ahora veremos el formato del archivo de salida que contiene las series sintéticas.
# Se carga el set de datos simulados
simulated_climate <- readr::read_csv(here::here('output_data/local/simulated_local_conditional.csv')) ##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## realization = col_double(),
## station_id = col_double(),
## point_id = col_double(),
## longitude = col_double(),
## latitude = col_double(),
## date = col_date(format = ""),
## tmax = col_double(),
## tmin = col_double(),
## prcp = col_double()
## )| realization | station_id | point_id | longitude | latitude | date | tmax | tmin | prcp |
|---|---|---|---|---|---|---|---|---|
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-01 | 28.05965 | 17.69135 | 7.534622 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-02 | 27.32783 | 16.20019 | 0.000000 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-03 | 35.01025 | 13.35766 | 0.000000 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-04 | 34.18990 | 18.04710 | 0.000000 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-05 | 30.59661 | 20.62692 | 0.000000 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-06 | 37.15162 | 17.44571 | 0.000000 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-07 | 35.65166 | 18.36680 | 0.000000 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-08 | 28.28271 | 13.77686 | 0.000000 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-09 | 29.66045 | 13.73328 | 0.000000 |
| 1 | 87448 | 1 | 5001636 | 6256577 | 2010-01-10 | 28.07379 | 14.17866 | 11.043809 |
Al utilizar totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima, las series generadas capturan los variaciones de baja frecuencia que se observan en la serie histórica. A continuación se muestra una Figura que ilustra lo anterior para la temperatura máxima.
ggplot2::ggplot() +
ggplot2::geom_boxplot(data = simulated_climate %>%
dplyr::mutate(month = lubridate::month(date),
year = lubridate::year(date),
date = as.Date(paste0(year, '-', month, '-', 15L))) ,
ggplot2::aes(x = date, y = tmax, group = date, fill = 'DarkOrange'), alpha = 0.1) +
ggplot2::geom_line(data = climate %>%
dplyr::filter(date > as.Date('2010-01-01')) %>%
dplyr::mutate(month = lubridate::month(date),
year = lubridate::year(date)) %>%
dplyr::group_by(year, month) %>%
dplyr::summarise(tmax = median(tmax, na.rm = TRUE)) %>%
dplyr::mutate(date = as.Date(paste0(year, '-', month, '-', 15L))),
ggplot2::aes(x = date, y = tmax, group = 1, color = 'DarkOrange')) +
ggplot2::theme_bw() +
ggplot2::labs(x = 'Fecha', y = 'Temperatura máxima [°C]') +
ggplot2::scale_x_date(breaks = '6 months', labels = scales::date_format("%m-%y")) +
ggplot2::scale_fill_discrete(name = "", labels = c("Simulado")) +
ggplot2::scale_color_discrete(name = "", labels = c("Observado")) +
ggplot2::theme(legend.position = 'bottom')## `summarise()` regrouping output by 'year' (override with `.groups` argument)
Las cajas corresponde a las distintas realizaciones agregadas a escala mensual y la línea naranja corresponde a la temperatura media mensual calculada para los años 2010 a 2019. Se observa como las cajas suben y bajan al ritmo de la media observada y como capturan los pequeños cambios que ocurren en un año específico pero no en otros.
Series sintéticas correlacionadas espacialmente
Una tercera alternativa para la generación de datos sobre estaciones meteorológicas combina las antes mostradas pero generando el tiempo local con métodos que contemplan la autocorrelación espacial. Para utilizar esta alternativa se deben disponer de más de una sola estación porque de otro modo no se podrían calcular los parámetros de los variogramas necesarios para la generación del tiempo local. Mientras más estaciones haya mejor, pero si pueden generar campos espaciales confiables con alrededor de 10 puntos.
Crear archivos de entrada
Para este ejemplo se utilizan datos de varias estaciones pero el formato de los mismos es igual a los mostrados anteriormente.
Los archivos necesarios son:
- stations.csv
- climate.csv
Los datos meteorológicos se dividen en dos archivos separados: stations.csv y climate.csv. Los nombres de los mismos no deben ser necesariamente iguales a los usados aquí.
Los metadatos de las estaciones se alojan en el archivo stations.csv. Este archivo contiene la información de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:
- station_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integer
- latitude: latitud en grados decimales. La variable debe ser de tipo double
- longitude: longitud en grados decimales. La variable debe ser de tipo double
La tabla puede tener más variables pero sólo se necesitan las anteriores.
A continuación se muestran la primera fila del dataset y los tipo de datos de cada una de las variables.
| station_id | nombre | lat_dec | lon_dec | elev |
|---|---|---|---|---|
| 86330 | Artigas | -30.4000 | -56.5100 | 120.38 |
| 86350 | Rivera | -30.9000 | -55.5400 | 241.94 |
| 86360 | Salto | -31.4300 | -57.9800 | 41.00 |
| 86430 | Paysandú | -32.3500 | -58.0400 | 61.12 |
| 86440 | Melo | -32.3700 | -54.1900 | 100.36 |
| 86460 | Paso de los Toros | -32.8000 | -56.5300 | 75.48 |
| 86490 | Mercedes | -33.2500 | -58.0700 | 17.01 |
| 86560 | Colonia | -34.4500 | -57.7700 | 18.00 |
| 86565 | Rocha | -34.4900 | -54.3100 | 18.16 |
| 86580 | Aeropuerto Carrasco | -34.8300 | -56.0100 | 32.88 |
| 90000001 | Las Brujas | -34.6719 | -56.3400 | 32.00 |
| 90000002 | La Estanzuela | -34.3372 | -57.6922 | 72.00 |
| 90000003 | Tacuarembó | -31.7089 | -55.8267 | 115.00 |
| 90000004 | Treinta y Tres | -33.2750 | -54.1722 | 22.00 |
| 90000005 | Salto Grande | -31.2728 | -57.8908 | 47.00 |
El objeto stations debe ser convertido de tibble a sf. El sistema de referencia espacial debe ser planar. No es necesario un sistema de referencia espacial en particular, solamente las coordenadas deben estar expresadas en metros.
# Convertimos el objeto stations a sf y se convierte su proyección de WGS 1984 a
# UTM Zona 21 S
stations %<>%
sf::st_as_sf(coords = c("lon_dec", "lat_dec"), crs = 4326) %>%
sf::st_transform(crs = 32721)En el siguiente mapa muestra la distribución de las estaciones meteorológicas usadas en el ejemplo.
# Convertir el objeto con las estaciones a coordenadas geográficas
stations_geo <- stations %>%
sf::st_transform(crs = 4326) %>%
dplyr::mutate(lat = sf::st_coordinates(.)[,2],
lon = sf::st_coordinates(.)[,1])
# Creación del mapa con la distribución de las estaciones
leaflet::leaflet(data = stations_geo) %>%
leaflet::addTiles() %>%
leaflet::setView(lat = -33, lng = -56, zoom = 5) %>%
leaflet::addCircleMarkers(lat = ~lat, lng = ~lon, radius = 3,
fillColor = "#377eb8",
color = "#377eb8",
stroke = FALSE, fillOpacity = 1, opacity = 1,
popup = ~sprintf("<b>%s (%d)</b><br>Lat.: %.3f<br>Lon.: %.3f<br>Elev: %.0f m",
nombre, station_id, lat, lon, elev)) La información climática se aloja en el archivo climate.csv. Este archivo contiene los datos de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:
date: fecha del dato. La variable debe ser de tipo datestation_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integerprcp: datos diarios de precipitación La variable debe ser de tipo doubletmax: datos diarios de temperatura máxima. La variable debe ser de tipo doubletmin: datos diarios de temperatura mínima. La variable debe ser de tipo double
A continuación se muestran las primeras cinco filas del dataset y los tipos de datos de cada una de las variables.
| date | station_id | tmax | tmin | prcp |
|---|---|---|---|---|
| 1981-01-01 | 86330 | 34.6 | 21.7 | 0.0 |
| 1981-01-02 | 86330 | 33.8 | 21.2 | 19.5 |
| 1981-01-03 | 86330 | 28.3 | 19.4 | 0.0 |
| 1981-01-04 | 86330 | 30.3 | 17.4 | 0.0 |
| 1981-01-05 | 86330 | 33.4 | 17.4 | 21.0 |
| 1981-01-06 | 86330 | 32.8 | 20.4 | 8.0 |
| 1981-01-07 | 86330 | 30.6 | 17.9 | 0.0 |
| 1981-01-08 | 86330 | 30.2 | 19.1 | 0.0 |
| 1981-01-09 | 86330 | 31.3 | 20.2 | 0.0 |
| 1981-01-10 | 86330 | 30.2 | 20.6 | 0.0 |
En total se utilizan 15 estaciones meteorologicas, 10 provistas por INUMET (Instituyo Uruguayo de Mteorología) y 5 por INIA (Instituto Nacional de Investigación Agropecuaria).
## [1] 86330 86350 86360 86430 86440 86460 86490 86560
## [9] 86565 86580 90000001 90000002 90000003 90000004 90000005Ajuste de los modelos
El ajuste de los modelos es igual que para los dos ejemplos antes mostrados. Se utiliza la función control_fit para la congifuración del ajuste y local_calibrate para realizar el ajuste.
control_fit <-gamwgen::local_fit_control(
prcp_occurrence_threshold = 0.1, # Umbral para la definición de días húmedos
avbl_cores = 6, # Cantidad de núcleos disponibles
planar_crs_in_metric_coords = 32721) # Sistema de referencia espacial (en metros)Este ejemplo se realizará utilizando covriables trimestrales pero es válido para series estacionarias.
# Agregación de valores diarios
seasonal_covariates <- gamwgen::summarise_seasonal_climate(climate, umbral_faltantes = 0.2)
# Se muestran las primeras cinco filas
knitr::kable(seasonal_covariates[1:10,])| station_id | year | season | seasonal_prcp | seasonal_tmax | seasonal_tmin |
|---|---|---|---|---|---|
| 86330 | 1981 | 1 | 458.9 | 30.32778 | 19.151111 |
| 86330 | 1981 | 2 | 370.4 | 26.44674 | 14.951087 |
| 86330 | 1981 | 3 | 211.2 | 19.76630 | 9.120652 |
| 86330 | 1981 | 4 | 272.0 | 24.73516 | 13.093407 |
| 86330 | 1982 | 1 | 465.9 | 30.33667 | 17.764444 |
| 86330 | 1982 | 2 | 229.0 | 26.74130 | 14.385870 |
| 86330 | 1982 | 3 | 365.4 | 19.48913 | 10.109783 |
| 86330 | 1982 | 4 | 651.7 | 24.68242 | 13.048352 |
| 86330 | 1983 | 1 | 609.0 | 31.39222 | 19.440000 |
| 86330 | 1983 | 2 | 471.2 | 23.80761 | 13.631522 |
En la siguiente Figura se muestra la variación de la precipitación acumulada estival para las 15 estaciones utilizadas.
ggplot2::ggplot(data = seasonal_covariates %>%
dplyr::filter(season == 1),
ggplot2::aes(x = year, y = seasonal_prcp)) +
ggplot2::geom_line() +
ggplot2::facet_wrap(~station_id, ncol = 3) +
ggplot2::theme_bw() +
ggplot2::labs(x = 'Años', y = 'Precipitación estival [mm]')
Se observa en la figura que si bien hay cierta coincidencia, no todos los picos y valles ocurren en el mismo momento. Utilizar un método que genere el tiempo local considerando estas variaciones espaciales es muy útil para representar feahacientemente los patrones espciales y temporales.
# Al correr la función se realiza el ajuste de los cuatro modelos para cada una de
# las estaciones. En este caso, por cuestiones de tiempo a cargar un objeto ya precalculado.
# Si el usuario desea correrlo deberá ver la nota anterior.
gamgen_fit <- gamwgen::local_calibrate(climate = climate, # Registro histórico de variables meteorológicas
stations = stations, # Estaciones meteorológicas
seasonal_covariates = seasonal_covariates, # Totales trimestrales de precipitación
control = control_fit, # Objeto de control
verbose = FALSE) # Impresión de mensajes en la consola.Para esta demostración, cargamos el objeto con el ajuste del modelo ya realizado.
# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('input_data/local/fit_local_correlated_weather.RData')) {
fs::file_copy(system.file('/autorun/local', "fit_local_correlated_weather.RData", package = "gamwgen"),
new_path = 'input_data/local/fit_local_correlated_weather.Rdata')
}
# Cargamos el archivo recientemente creado
load('input_data/local/fit_local_correlated_weather.RData')
# Clase del objeto con el ajuste del generador
class(gamgen_fit)## [1] "gamwgen"## [1] "control" "stations" "climate"
## [4] "seasonal_covariates" "crs_used_to_fit" "start_climatology"
## [7] "fitted_models" "models_data" "models_residuals"
## [10] "statistics_threshold" "exec_times"Dentro del objeto se guardan todo lo necesario para la simulación así como información accesoria.
control: copia de la configuración usada para calibrar el generadorstations: estaciones meteorológicas utilizadas para la calibraciónclimate: datos climáticos de cada uno de las estacionesseasonal_covariates: series temporales de totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima.crs_used_to_fit: sistema de referencia espacial usado para proyectarstart_climatology: climatología diaria de cada una de las variables de entrada.fitted_models: modelos ajustados, uno para cada variable: temperaturas máxima y mínima y ocurrencia y montos de precipitación.models_data: datos usados efectivamente usados para ajustar los modelos (sin NAs)models_residuals: residuos de cada uno de los modelos. Es decir, la diferencia entre el valor ajustado por el modelo (clima local) y el valor observado en el día istatistics_threshold: umbrales de amplitud térmica diaria por mes. Si la amplitud simulada está fuera de este rango, se repetirá la simulación para ese día a los fines de mantener la consistencia entre variablesexec_times: tiempo de ejecución de cada una de las etapas del ajuste Al inspeccionar el objeto con los resultados se puede ver que cada estación meteorológica tiene su propia sublista donde se almacenan los modelos para cada una de ellas.
## [1] "86330" "86350" "86360" "86430" "86440" "86460"
## [7] "86490" "86560" "86565" "86580" "90000001" "90000002"
## [13] "90000003" "90000004" "90000005"Cada uno de los GAMs ajustados se almacenan en el objeto gamgen_fit y pueden ser evaluados con la función summary().
##
## Family: gaussian
## Link function: identity
##
## Formula:
## tmax ~ s(tmax_prev, tmin_prev, k = 50) + s(prcp_occ, bs = "re") +
## s(prcp_occ_prev, bs = "re") + s(doy, bs = "cc", k = 30) +
## s(SX1, SN1, k = 20) + s(SX2, SN2, k = 20) + s(SX3, SN3, k = 20) +
## s(SX4, SN4, k = 20)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.18758 0.03397 770.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(tmax_prev,tmin_prev) 32.4286 40.704 235.57 <2e-16 ***
## s(prcp_occ) 0.9979 1.000 520.81 <2e-16 ***
## s(prcp_occ_prev) 0.9987 1.000 866.16 <2e-16 ***
## s(doy) 10.3038 28.000 33.61 <2e-16 ***
## s(SX1,SN1) 3.6312 4.458 26.92 <2e-16 ***
## s(SX2,SN2) 2.0001 2.000 36.90 <2e-16 ***
## s(SX3,SN3) 4.3491 5.631 21.09 <2e-16 ***
## s(SX4,SN4) 2.0007 2.001 40.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.778 Deviance explained = 77.9%
## fREML = 35496 Scale est. = 8.8219 n = 14128La función summary permite analizar los resultados del ajuste de cada uno de los modelos. Para el caso del modelo de temperatura máxima podemos ver la fórmula del GAM en la parte superior bajo el apartado Formula y la significancia de cada uno de los términos del modelo en la tabla inmediatamente inferior. En este configuracion, al incluirle variables estacionales, se incoporan nuevos términos el modelo como SX1, SX2, SX3, SX4 y SN1, SN2, SN3, SN4, que corresponden a variables dummy de las medias trimestrales de temperatura máxima y mínima, respectivamente. Se puede observar que todos los términos son altamente significativos. También se incluyen como pruebas de bondad del ajuste el porcentaje de la varianza explicada por el modelo y el valor de R-ajustado.
El paquete gratia (Simpson, Gavin (2018)) es muy útil para la visualización de los ajustes de manera gráfica. Esta función toma los diagnósticos por defecto de la función gam.check del paquete mgcv (Wood, Simon (2001)) y los convierte en gráficos de la librería ggplot2 (Wickham, Hadley(2011)) que son visualmente muy atractivos. La figura está dividida en cuatro paneles, cada uno mostrando un diagnóstico diferente. En el panel superior izquierdo se muestra un Q-Q Plot de los residuos. Si el modelo ha ajustado satisfactoriamente los datos, los puntos deberían estar sobre la recta 1:1 demarcada en rojo. El panel superior derecho muestra los residuos vs el predictor lineal. El objetivo de este diagnóstico es que los residuos se distribuyan al azar y que no tengan un patrón claro. La presencia de patrones en los residuos indicaría que el modelo no ha explicado algun componente importante de la variabilidad de los datos. En el panel inferior izquierdo se muestra un histograma de los residuos siendo deseable que lo mismos tengan una distribución proxima a la Normal. El último gráfico en el panel inferior derecho muestra una gráfico de dispersión entre los valores ajustados y observados.
Estos diagnóstico exploratorios puede aplicarse a cada uno de los modelos ajustados por la función local_calibrate.
Con la creación del objeto gamgen_fit finaliza el proceso de ajuste y ya se pueden generar series sintéticas para la o las estaciones usadas para calibrar el generador.
Generación de series
La generación de series sigue la misma estructura anterior, una función para configurar la generación y otra que realiza la generación propiamente dicha.
Los argumentos de la función de control son:
nsim: cantidad de simulaciones a realizar. Se debe ingresar un valor entero mayor o igual a 1.seed: semilla. Se debe ingresar cualquier numero entero. No es necesario recordarlo porque se guarda junto a los resultados.avbl_cores: cantidad de núcleos disponibles para la paralelización.use_spatially_correlated_noise: utilizar la generación estocástica espacialmente correlacionada. Esta opción sólo es válida si en el ajuste y en la simulación se usaron más de cinco estaciones meteorológicas diferentes. Con un menor número no es posible calcular los variogramas necesarios para la generación de los campos aleatorios. Se debe introducir un boolean (TRUE or FALSE).use_temporary_files_to_save_ram: si se simulan muchas realizaciones o los recursos informáticos son escasos, esta opción permite guardar los resultados de cada una de las realizaciones en el disco liberando memoria RAM que quedará disponible para generar nuevas simulaciones. Al finalizar la generación todos los archivos se combinan en uno único. Se debe introducir un boolean (TRUE or FALSE).use_temporary_files_to_save_ram: esta opción permite eliminar los archivos temporales creados para ahorrar RAM luego de terminar la generación de todas las simulaciones. Se debe introducir un boolean (TRUE or FALSE).
control_sim <- gamwgen::local_simulation_control(
nsim = 10, # Cantidad de simulaciones a realizar
seed = 1234, # Semilla para que los resultados sean reproducibles
avbl_cores = 6, # Cantidad de núcleos disponibles a utilizar
use_spatially_correlated_noise = TRUE, # Usar modelo de ruido espacialmente correlacionado
use_temporary_files_to_save_ram = TRUE, # Guardar resultados intermedios para ahorrar RAM
remove_temp_files_used_to_save_ram = TRUE) # Borrar los resultados intermedios creados anteriormenteLuego se procede a la simulación de datos meteorológicos. Los argumentos de la función de simulación son:
model: objeto con el resultado de la funciónlocal_calibrate()simulation_locations: objeto tiposfcon la ubicación de las estaciones a simular. Las estaciones usadas deben haber sido incluidas en el proceso de ajuste. No es necesario que todas estén presentes, se pueden generar series solo sobre algunas de ellas.start_date: fecha de comienzo de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de comienzo es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, ni tampoco posterior. Se debe introducir una fecha en formato dateend_date: fecha de fin de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de fin es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, tampoco puede ser posterior. Se debe introducir una fecha en formato datecontrol: objeto de control creado con la funcióncontrol_sim().output_folder: ruta al directorio donde se guardarán los resultados, tanto finales como intermedios.output_filename: nombre del archivo de salida. Para facilitar la interoperabilidad, el archivo generado es un archivo de texto en formato separado por comas (.csv)seasonal_covariates: datos agregados trimestrales. Si el ajuste se realizó con covariables, la generación también debe realizarse con ellas. Caso contrario se producirá un error. Se debe introducir un data frame con los valores agregados para las tres variables (precipitación y temperaturas máxima y mínima) pero no necesariamente deben ser los mismos a los utilizados en el ajuste. Si se desean simular tendencias de algún tipo, ya sea de un modelo de cambio climático o arbitrarias, se deben perturbar estas variables trimestrales e introducirlas aquí. Estas series si deben tener la misma longitud que el período a generarverbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza la generación de series para cada una de las estaciones.
# En este caso, por cuestiones de tiempo, vamos a cargar un objeto con los resultados de la simulación
simulated_climate <- gamwgen::local_simulation(model = gamgen_fit, # Objeto con los resultados del ajuste
simulation_locations = stations, # Estaciones para las cuales simular
start_date = as.Date('2019-01-01'), # Fecha de comienzo de las simulaciones
end_date = as.Date('2019-12-31'), # Fecha de fin de las simulaciones
control = control_sim, # Objeto con la configuración
output_folder = getwd(), # Directorio donde se guardarán los resultados
output_filename = 'simulations.csv', # Nombre del archivo de salida
seasonal_covariates = seasonal_covariates, # Covariables estacionales
verbose = FALSE) # Impresión de mensajes en la consolaEsta función produce dos tipos de resultados: una lista que permanece en el ambiente de R y los datos generados que son guardados como .csv en el directorio indicado precedentemente.
# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('output_data/local/simulatied_local_correlated_weather.RData')) {
fs::file_copy(system.file('/autorun/local', "simulatied_local_correlated_weather.RData", package = "gamwgen"),
new_path = 'output_data/local/simulatied_local_correlated_weather.RData')
}
# Cargamos el archivo recientemente creado
load('output_data/local/simulated_local_unconditional.RData')
# Clase del objeto con el ajuste del generador
class(simulated_climate)## [1] "list" "gamwgen.climate"## [1] "nsim"
## [2] "seed"
## [3] "realizations_seeds"
## [4] "simulation_points"
## [5] "output_file_with_results"
## [6] "output_file_fomart"
## [7] "rdata_file_with_fitted_stations_and_climate"
## [8] "exec_times"La lista contiene los siguientes objetos:
nsim: cantidad de realizaciones.seed: semilla general para toda la generación. Corresponde a la que se incluye en la función de control.realization_seeds: semillas para cada una de las realizaciones. Esto permite replicar los resultados.simulation_points: puntos donde se generaron las series sintéticas.output_fil_with_results: nombre del archivo con los resultados.output_file_format: tipo de archivo de salida, en este caso .csv.rdata_file_with_fitted_stations_and_climate: archivo.RDatacon los datos meteorológicos observados que fueron utilizados en el ajuste. También se incluyen los metadatos de cada uno de esos puntos.exec_times: tiempo de ejecución del la generación.
Ahora veremos el formato del archivo de salida que contiene las series sintéticas.
Para desmenuzar la generación del tiempo local se pueden correr algunas de las funciones que forman parte de local_simulation. Por ejemplo, del objeto gamgen_fit se extraen los residuos y con la función gamwgen:::generate_residuals_statistics se calculan los parámetros.
En la tabla anterior se muestran los parámetros necesarios para generar el tiempo local de las temperaturas máxima y mínima.
gen_noise_params <- gamwgen:::generate_month_params(
residuals = gamgen_fit$models_residuals,
observed_climate = gamgen_fit$models_data,
stations = gamgen_fit$stations)
# Ejemplo de variograma de residuos de temperatura máxima para días secos
gen_noise_params[[1]]$variogram_parameters$tmax_dry ## [1] 0.000000 6.978337 495.779099Para cada una de las variables, precipitación y temperaturas máxima y mínima, se ajusta un variograma por máxima verosimilitud que permite representar la variabilidad espacial de los residuos de los modelos GAMs, que corresponden al tiempo local. En el variograma se muestran los tres parámetros nugget, psill y range. El rango de los datos es de 500 km, lo que es lógico si se considera que los puntos están distribuidamente en la República Oriental del Uruguay.
En la siguiente FIgura se muestra la variabilidad espacial del tiempo local de temperatura del día 1 de enero de 2019 para días secos. Cabe mencionar que se crean también los mismos datos para los días húmedos y en función de la ocurrencia de precipitación se selecciona uno u otro.
RandomFields::RFoptions(printlevel = 0, spConform = FALSE)
temp_local <- control_sim$temperature_noise_generating_function(simulation_points = stations %>%
dplyr::mutate(latitude = sf::st_coordinates(.)[,2],
longitude = sf::st_coordinates(.)[,1]),
gen_noise_params = gen_noise_params,
month_number = 1L,
selector = 'Dry',
seed = NULL) %>%
dplyr::mutate(date = as.Date(paste0('2019-', '01', '-', '01'))) %>%
tidyr::gather(variable, valor, -c('date', 'geometry'))
# Mapa de Uruguay
uruguay <- raster::getData('GADM', country='URY', level=1) %>%
sf::st_as_sf(.) %>%
sf::st_transform(crs = 32721)
# Paleta de colores
colr <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(9, 'YlOrRd'))
# Quiebres de la escala de colores
quiebres <- c(-8.0, -6.0, -4.0,-2.0, 0, 2.0, 4.0, 6.0, 8.0)
# Etiquetas para los colores
etiquetas <- c('-8', '-6', '-4', '-2', '0', '2', '4', '6', '8')
# Etiquetas de las facetas
labs <- c("Temp. Máx", "Temp. Min.")
names(labs) <- c("tmax_residuals", "tmin_residuals")
ggplot2::ggplot(data = temp_local) +
ggplot2::geom_sf(ggplot2::aes(color = valor), size = 4) +
ggplot2::scale_color_gradientn(name = "Tiempo local [°C]",
colours = colr(9),
breaks = quiebres,
labels = etiquetas,
limits = c(-8, 8)) +
ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
ggplot2::facet_wrap(~variable, labeller = ggplot2::labeller(variable = labs)) +
ggplot2::theme_bw()
Se puede observar como el tiempo local tiene una estructura espacial que es representada a través del variograma y permite que los resultados para las distintas estaciones sean espacialmente consistentes. Esto quiere decir que estaciones meteorológicas cercanas tenderán a tener valores de tiempo local similares.
Generación de series sintéticas para una grilla regular
La generación de series sobre grillas regular es muy similar a lo mostrado anterioemente. Los fundamentos son los mismos, sólo los modelos estadísticos que modelan el clima local incoporan un nuevo término que permite incluir la dimensión espacial.
Crear archivos de entrada
Para este ejemplo se utilizan datos de varias estaciones pero el formato de los mismos es igual a los mostrados anteriormente.
Los archivos necesarios son:
- stations.csv
- climate.csv
Los datos meteorológicos se dividen en dos archivos separados: stations.csv y climate.csv. Los nombres de los mismos no deben ser necesariamente iguales a los usados aquí.
Los metadatos de las estaciones se alojan en el archivo stations.csv. Este archivo contiene la información de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:
- station_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integer
- latitude: latitud en grados decimales. La variable debe ser de tipo double
- longitude: longitud en grados decimales. La variable debe ser de tipo double
La tabla puede tener más variables pero sólo se necesitan las anteriores.
A continuación se muestran la primera fila del dataset y los tipo de datos de cada una de las variables.
| station_id | nombre | lat_dec | lon_dec | elev |
|---|---|---|---|---|
| 86330 | Artigas | -30.4000 | -56.5100 | 120.38 |
| 86350 | Rivera | -30.9000 | -55.5400 | 241.94 |
| 86360 | Salto | -31.4300 | -57.9800 | 41.00 |
| 86430 | Paysandú | -32.3500 | -58.0400 | 61.12 |
| 86440 | Melo | -32.3700 | -54.1900 | 100.36 |
| 86460 | Paso de los Toros | -32.8000 | -56.5300 | 75.48 |
| 86490 | Mercedes | -33.2500 | -58.0700 | 17.01 |
| 86560 | Colonia | -34.4500 | -57.7700 | 18.00 |
| 86565 | Rocha | -34.4900 | -54.3100 | 18.16 |
| 86580 | Aeropuerto Carrasco | -34.8300 | -56.0100 | 32.88 |
| 90000001 | Las Brujas | -34.6719 | -56.3400 | 32.00 |
| 90000002 | La Estanzuela | -34.3372 | -57.6922 | 72.00 |
| 90000003 | Tacuarembó | -31.7089 | -55.8267 | 115.00 |
| 90000004 | Treinta y Tres | -33.2750 | -54.1722 | 22.00 |
| 90000005 | Salto Grande | -31.2728 | -57.8908 | 47.00 |
El objeto stations debe ser convertido de tibble a sf. El sistema de referencia espacial debe ser planar. No es necesario un sistema de referencia espacial en particular, solamente las coordenadas deben estar expresadas en metros.
# Convertimos el objeto stations a sf y se convierte su proyección de WGS 1984 a
# UTM Zona 21 S
stations %<>%
sf::st_as_sf(coords = c("lon_dec", "lat_dec"), crs = 4326) %>%
sf::st_transform(crs = 32721)En el siguiente mapa muestra la distribución de las estaciones meteorológicas usadas en el ejemplo.
# Convertir el objeto con las estaciones a coordenadas geográficas
stations_geo <- stations %>%
sf::st_transform(crs = 4326) %>%
dplyr::mutate(lat = sf::st_coordinates(.)[,2],
lon = sf::st_coordinates(.)[,1])
# Creación del mapa con la distribución de las estaciones
leaflet::leaflet(data = stations_geo) %>%
leaflet::addTiles() %>%
leaflet::setView(lat = -33, lng = -56, zoom = 5) %>%
leaflet::addCircleMarkers(lat = ~lat, lng = ~lon, radius = 3,
fillColor = "#377eb8",
color = "#377eb8",
stroke = FALSE, fillOpacity = 1, opacity = 1,
popup = ~sprintf("<b>%s (%d)</b><br>Lat.: %.3f<br>Lon.: %.3f<br>Elev: %.0f m",
nombre, station_id, lat, lon, elev)) La información climática se aloja en el archivo climate.csv. Este archivo contiene los datos de las estaciones meteorológicas que serán usadas en el ajuste del modelo. Las variables que deben ser incluidas en la tabla son:
date: fecha del dato. La variable debe ser de tipo datestation_id: número unívoco para cada estación meteorológica. La variable debe ser de tipo integerprcp: datos diarios de precipitación La variable debe ser de tipo doubletmax: datos diarios de temperatura máxima. La variable debe ser de tipo doubletmin: datos diarios de temperatura mínima. La variable debe ser de tipo double
A continuación se muestran las primeras cinco filas del dataset y los tipos de datos de cada una de las variables.
| date | station_id | tmax | tmin | prcp |
|---|---|---|---|---|
| 1981-01-01 | 86330 | 34.6 | 21.7 | 0.0 |
| 1981-01-02 | 86330 | 33.8 | 21.2 | 19.5 |
| 1981-01-03 | 86330 | 28.3 | 19.4 | 0.0 |
| 1981-01-04 | 86330 | 30.3 | 17.4 | 0.0 |
| 1981-01-05 | 86330 | 33.4 | 17.4 | 21.0 |
| 1981-01-06 | 86330 | 32.8 | 20.4 | 8.0 |
| 1981-01-07 | 86330 | 30.6 | 17.9 | 0.0 |
| 1981-01-08 | 86330 | 30.2 | 19.1 | 0.0 |
| 1981-01-09 | 86330 | 31.3 | 20.2 | 0.0 |
| 1981-01-10 | 86330 | 30.2 | 20.6 | 0.0 |
En total se utilizan 15 estaciones meteorologicas, 10 provistas por INUMET (Instituyo Uruguayo de Meteorología) y 5 por INIA (Instituto Nacional de Investigación Agropecuaria).
## [1] 86330 86350 86360 86430 86440 86460 86490 86560
## [9] 86565 86580 90000001 90000002 90000003 90000004 90000005Ajuste de los modelos
El ajuste de los modelos es igual que para los dos ejemplos antes mostrados. Se utiliza la función control_fit para la congifuración del ajuste y spatial_calibrate para realizar el ajuste.
control_fit <-gamwgen::spatial_fit_control(
prcp_occurrence_threshold = 0.1, # Umbral para la definición de días húmedos
avbl_cores = 6, # Cantidad de núcleos disponibles
planar_crs_in_metric_coords = 32721) # Sistema de referencia espacial (en metros)Este ejemplo se realizará utilizando covriables trimestrales pero es válido para series estacionarias.
# Agregación de valores diarios
seasonal_covariates <- gamwgen::summarise_seasonal_climate(climate, umbral_faltantes = 0.2)
# Se muestran las primeras cinco filas
knitr::kable(seasonal_covariates[1:10,])| station_id | year | season | seasonal_prcp | seasonal_tmax | seasonal_tmin |
|---|---|---|---|---|---|
| 86330 | 1981 | 1 | 458.9 | 30.32778 | 19.151111 |
| 86330 | 1981 | 2 | 370.4 | 26.44674 | 14.951087 |
| 86330 | 1981 | 3 | 211.2 | 19.76630 | 9.120652 |
| 86330 | 1981 | 4 | 272.0 | 24.73516 | 13.093407 |
| 86330 | 1982 | 1 | 465.9 | 30.33667 | 17.764444 |
| 86330 | 1982 | 2 | 229.0 | 26.74130 | 14.385870 |
| 86330 | 1982 | 3 | 365.4 | 19.48913 | 10.109783 |
| 86330 | 1982 | 4 | 651.7 | 24.68242 | 13.048352 |
| 86330 | 1983 | 1 | 609.0 | 31.39222 | 19.440000 |
| 86330 | 1983 | 2 | 471.2 | 23.80761 | 13.631522 |
# Al correr la función se realiza el ajuste de los cuatro modelos para cada una de
# las estaciones. En este caso, por cuestiones de tiempo a cargar un objeto ya precalculado.
# Si el usuario desea correrlo deberá ver la nota anterior.
gamgen_fit <- gamwgen::spatial_calibrate(climate = climate, # Registro histórico de variables meteorológicas
stations = stations, # Estaciones meteorológicas
seasonal_covariates = seasonal_covariates, # Totales trimestrales de precipitación
control = control_fit, # Objeto de control
verbose = FALSE) # Impresión de mensajes en la consola.Para esta demostración, cargamos el objeto con el ajuste del modelo ya realizado.
# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('input_data/spatial/fit_spatial_conditional.RData')) {
fs::file_copy(system.file('/autorun/spatial', "fit_spatial_conditional.RData", package = "gamwgen"),
new_path = 'input_data/spatial/fit_spatial_conditional.Rdata')
}
# Cargamos el archivo recientemente creado
load('input_data/spatial/fit_spatial_conditional.RData')
# Clase del objeto con el ajuste del generador
class(gamgen_fit)## [1] "gamwgen"## [1] "control" "stations" "climate"
## [4] "seasonal_covariates" "crs_used_to_fit" "start_climatology"
## [7] "fitted_models" "models_data" "models_residuals"
## [10] "statistics_threshold" "exec_times"Dentro del objeto se guardan todo lo necesario para la simulación así como información accesoria.
control: copia de la configuración usada para calibrar el generadorstations: estaciones meteorológicas utilizadas para la calibraciónclimate: datos climáticos de cada uno de las estacionesseasonal_covariates: series temporales de totales trimestrales de precipitación y medias trimestrales de temperaturas máxima y mínima.crs_used_to_fit: sistema de referencia espacial usado para proyectarstart_climatology: climatología diaria de cada una de las variables de entrada.fitted_models: modelos ajustados, uno para cada variable: temperaturas máxima y mínima y ocurrencia y montos de precipitación.models_data: datos usados efectivamente usados para ajustar los modelos (sin NAs)models_residuals: residuos de cada uno de los modelos. Es decir, la diferencia entre el valor ajustado por el modelo (clima local) y el valor observado en el día istatistics_threshold: umbrales de amplitud térmica diaria por mes. Si la amplitud simulada está fuera de este rango, se repetirá la simulación para ese día a los fines de mantener la consistencia entre variablesexec_times: tiempo de ejecución de cada una de las etapas del ajuste
A diferencia de la configuración que ajusta distintos modelos para cada estación meteorológica, en te caso se ajusta un sólo modelos para toda la región de interés. Al visualizar lo que está almacenado en la sublista fitted_models se ve que solo hay cuatro GAMs, dos para la precipitación y dos para las temperaturas máxima y mínima.
## [1] "prcp_occ_fit" "prcp_amt_fit" "tmax_fit" "tmin_fit"Cada uno de los GAMs ajustados se almacenan en el objeto gamgen_fit y pueden ser evaluados con la función summary().
##
## Family: gaussian
## Link function: identity
##
## Formula:
## tmax ~ te(tmax_prev, tmin_prev, longitude, latitude, d = c(2,
## 2), k = length(unique_stations)) + te(type_day, longitude,
## latitude, d = c(1, 2), bs = c("re", "tp"), k = length(unique_stations)) +
## te(type_day_prev, longitude, latitude, d = c(1, 2), bs = c("re",
## "tp"), k = length(unique_stations)) + te(doy, longitude,
## latitude, d = c(1, 2), bs = c("cc", "tp"), k = length(unique_stations)) +
## te(SX1, SN1, longitude, latitude, d = c(2, 2), k = length(unique_stations)) +
## te(SX2, SN2, longitude, latitude, d = c(2, 2), k = length(unique_stations)) +
## te(SX3, SN3, longitude, latitude, d = c(2, 2), k = length(unique_stations)) +
## te(SX4, SN4, longitude, latitude, d = c(2, 2), k = length(unique_stations))
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.91 29.52 0.878 0.38
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## te(tmax_prev,tmin_prev,longitude,latitude) 76.902 86.484 1267.40 <2e-16 ***
## te(type_day,longitude,latitude) 27.680 29.000 161.22 <2e-16 ***
## te(type_day_prev,longitude,latitude) 14.703 15.000 726.24 <2e-16 ***
## te(doy,longitude,latitude) 47.907 195.000 77.44 <2e-16 ***
## te(SX1,SN1,longitude,latitude) 18.355 23.456 75.14 <2e-16 ***
## te(SX2,SN2,longitude,latitude) 6.005 6.009 184.32 <2e-16 ***
## te(SX3,SN3,longitude,latitude) 30.642 36.395 48.54 <2e-16 ***
## te(SX4,SN4,longitude,latitude) 6.002 6.004 193.99 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.779 Deviance explained = 77.9%
## fREML = 5.221e+05 Scale est. = 8.9846 n = 207269La función summary permite analizar los resultados del ajuste de cada uno de los modelos. Para el caso del modelo de temperatura máxima podemos ver la fórmula del GAM en la parte superior bajo el apartado Formula y la significancia de cada uno de los términos del modelo en la tabla inmediatamente inferior. En este configuracion, al incluirle variables estacionales, se incoporan nuevos términos el modelo como SX1, SX2, SX3, SX4 y SN1, SN2, SN3, SN4, que corresponden a variables dummy de las medias trimestrales de temperatura máxima y mínima, respectivamente. Se puede observar que todos los términos son altamente significativos. También se incluyen como pruebas de bondad del ajuste el porcentaje de la varianza explicada por el modelo y el valor de R-ajustado.
El paquete gratia (Simpson, Gavin (2018)) es muy útil para la visualización de los ajustes de manera gráfica. Esta función toma los diagnósticos por defecto de la función gam.check del paquete mgcv (Wood, Simon (2001)) y los convierte en gráficos de la librería ggplot2 (Wickham, Hadley(2011)) que son visualmente muy atractivos. La figura está dividida en cuatro paneles, cada uno mostrando un diagnóstico diferente. En el panel superior izquierdo se muestra un Q-Q Plot de los residuos. Si el modelo ha ajustado satisfactoriamente los datos, los puntos deberían estar sobre la recta 1:1 demarcada en rojo. El panel superior derecho muestra los residuos vs el predictor lineal. El objetivo de este diagnóstico es que los residuos se distribuyan al azar y que no tengan un patrón claro. La presencia de patrones en los residuos indicaría que el modelo no ha explicado algun componente importante de la variabilidad de los datos. En el panel inferior izquierdo se muestra un histograma de los residuos siendo deseable que lo mismos tengan una distribución proxima a la Normal. El último gráfico en el panel inferior derecho muestra una gráfico de dispersión entre los valores ajustados y observados.
Estos diagnóstico exploratorios puede aplicarse a cada uno de los modelos ajustados por la función spatial_calibrate.
Con la creación del objeto gamgen_fit finaliza el proceso de ajuste y ya se pueden generar series sintéticas para la o las estaciones usadas para calibrar el generador.
Generación de series
La generación de series sigue la misma estructura anterior, una función para configurar la generación y otra que realiza la generación propiamente dicha.
Los argumentos de la función de control son:
nsim: cantidad de simulaciones a realizar. Se debe ingresar un valor entero mayor o igual a 1.seed: semilla. Se debe ingresar cualquier numero entero. No es necesario recordarlo porque se guarda junto a los resultados.avbl_cores: cantidad de núcleos disponibles para la paralelización.use_spatially_correlated_noise: utilizar la generación estocástica espacialmente correlacionada. Esta opción sólo es válida si en el ajuste y en la simulación se usaron más de cinco estaciones meteorológicas diferentes. Con un menor número no es posible calcular los variogramas necesarios para la generación de los campos aleatorios. Se debe introducir un boolean (TRUE or FALSE).use_temporary_files_to_save_ram: si se simulan muchas realizaciones o los recursos informáticos son escasos, esta opción permite guardar los resultados de cada una de las realizaciones en el disco liberando memoria RAM que quedará disponible para generar nuevas simulaciones. Al finalizar la generación todos los archivos se combinan en uno único. Se debe introducir un boolean (TRUE or FALSE).use_temporary_files_to_save_ram: esta opción permite eliminar los archivos temporales creados para ahorrar RAM luego de terminar la generación de todas las simulaciones. Se debe introducir un boolean (TRUE or FALSE).
control_sim <- gamwgen::spatial_simulation_control(
nsim = 10, # Cantidad de simulaciones a realizar
seed = 1234, # Semilla para que los resultados sean reproducibles
avbl_cores = 6, # Cantidad de núcleos disponibles a utilizar
use_spatially_correlated_noise = TRUE, # Usar modelo de ruido espacialmente correlacionado
use_temporary_files_to_save_ram = TRUE, # Guardar resultados intermedios para ahorrar RAM
remove_temp_files_used_to_save_ram = TRUE) # Borrar los resultados intermedios creados anteriormente# Agregación de valores diarios
seasonal_covariates <- gamwgen::summarise_seasonal_climate(climate, umbral_faltantes = 0.2)# Shapefile de Uruguay
uruguay <- raster::getData('GADM', country='URY', level=1) %>%
sf::st_as_sf(.) %>%
sf::st_transform(crs = 32721)
# Creación de la grilla sobre la que se simulará
grilla_simulacion_centers <- uruguay %>%
dplyr::select(geometry) %>%
sf::st_transform(gamgen_fit$crs_used_to_fit) %>%
sf::st_make_grid(cellsize = 25000, what = 'centers') %>%
sf::st_sf(grid_id = 1:length(.), geometry = .) %>%
sf::st_intersection(uruguay)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesggplot2::ggplot() +
ggplot2::geom_sf(data = grilla_simulacion_centers) +
ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
ggplot2::theme_bw()
Al utilizar covariables trimestrales es necesario que éstas se interpolen para cada uno de los puntos de la grilla. A continuación se muestran dos ejemplos, uno para precipitación acumulada para el trimestre estival de 2019 y la temperatura máxima media del mismo período.
simulation_dates <-
tibble::tibble(date = seq.Date(from = as.Date('2019-01-01'),
to = as.Date('2019-01-31'),
by = "days")) %>%
dplyr::mutate(year = lubridate::year(date),
month = lubridate::month(date),
season = lubridate::quarter(date, fiscal_start = 12))
grilla_simulacion_centers <- uruguay %>%
sf::st_transform(gamgen_fit$crs_used_to_fit) %>%
sf::st_make_grid(cellsize = 25000, what = 'centers') %>%
sf::st_sf(grid_id = 1:length(.), geometry = .) %>%
sf::st_intersection(uruguay) %>%
dplyr::select(geometry) %>%
dplyr::mutate(latitude = sf::st_coordinates(.)[,2],
longitude = sf::st_coordinates(.)[,1])## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesseasonal_covariates <-
gamwgen:::get_covariates(model = gamgen_fit,
simulation_points = grilla_simulacion_centers,
seasonal_covariates,
simulation_dates,
control = control_sim)## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output## Warning in proj4string(d$data): CRS object has comment, which is lost in output## Warning in proj4string(newdata): CRS object has comment, which is lost in output
## Warning in proj4string(newdata): CRS object has comment, which is lost in output## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output## Warning in proj4string(d$data): CRS object has comment, which is lost in output## Warning in proj4string(newdata): CRS object has comment, which is lost in output
## Warning in proj4string(newdata): CRS object has comment, which is lost in output## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output## Warning in proj4string(input_data): CRS object has comment, which is lost in
## output## Warning in proj4string(new_data): CRS object has comment, which is lost in
## output## Warning in proj4string(d$data): CRS object has comment, which is lost in output## Warning in proj4string(newdata): CRS object has comment, which is lost in output
## Warning in proj4string(newdata): CRS object has comment, which is lost in outputaux <- gamwgen:::sf2raster(seasonal_covariates %>%
dplyr::filter(year == 2019, season == 1),
variable = 'seasonal_prcp') %>%
raster::as.data.frame(., xy = TRUE) %>%
dplyr::rename(., 'longitude' = 'x', 'latitude' = 'y') %>%
dplyr::mutate(., x = longitude, y = latitude) %>%
sf::st_as_sf(., coords = c( 'x', 'y')) %>%
sf::st_set_crs(., 32721) %>%
sf::st_intersection(uruguay) %>%
dplyr::select(latitude, longitude, z, geometry) ## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriescolores <- colorRampPalette(c("#cccccc", "#7bccc4", "#2b8cbe", "#0868ac",
"#084081", "#807dba", "#6a51a3", "#54278f", "#3f007d", "#e7298a",
"#ce1256", "#67001f"))
quiebres <- c(0.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0)
etiquetas <- c('0', '5', '10', '20', '30', '40', '50', '60', '70', '80', '90', '100', '110', '120', '130', '140', '150', '160')
ggplot(data = aux,
ggplot2::aes(x = longitude, y = latitude, fill = z)) +
geom_tile(colour="black") +
scale_fill_gradientn(name = "Precipitation [mm]",
colours = colores(17),
breaks = quiebres,
labels = etiquetas) +
ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
ggplot2::theme_bw() +
ggplot2::theme(panel.background = element_rect(fill = "ivory"),
axis.text = element_text(size = 11, colour = 1),
legend.key.height = unit(2.5, "cm"))
aux <- gamwgen:::sf2raster(seasonal_covariates %>%
dplyr::filter(year == 2019, season == 1),
variable = 'seasonal_tmax') %>%
raster::as.data.frame(., xy = TRUE) %>%
dplyr::rename(., 'longitude' = 'x', 'latitude' = 'y') %>%
dplyr::mutate(., x = longitude, y = latitude) %>%
sf::st_as_sf(., coords = c( 'x', 'y')) %>%
sf::st_set_crs(., 32721) %>%
sf::st_intersection(uruguay) %>%
dplyr::select(latitude, longitude, z, geometry) ## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriescolr <- colorRampPalette(RColorBrewer::brewer.pal(9, 'YlOrRd'))
quiebres <- c(0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0)
etiquetas <- c('0', '5', '10', '15', '20', '25', '30', '35', '40', '45')
ggplot(data = aux,
ggplot2::aes(x = longitude, y = latitude, fill = z)) +
geom_tile(colour="black") +
scale_fill_gradientn(name = "Temperatura máxima [°C]",
colours = colr(9),
breaks = quiebres,
labels = etiquetas) +
ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
ggplot2::theme_bw() +
ggplot2:: labs(title = "", x = "", y = "") +
ggplot2::theme(panel.background = element_rect(fill = "ivory"),
axis.text = element_text(size = 11, colour = 1))
Luego se procede a la simulación de datos meteorológicos. Los argumentos de la función de simulación son:
model: objeto con el resultado de la funciónlocal_calibrate()simulation_locations: objeto tiposfcon la ubicación de las estaciones a simular. Las estaciones usadas deben haber sido incluidas en el proceso de ajuste. No es necesario que todas estén presentes, se pueden generar series solo sobre algunas de ellas.start_date: fecha de comienzo de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de comienzo es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, ni tampoco posterior. Se debe introducir una fecha en formato dateend_date: fecha de fin de la generación de series sintéticas. Si no se incluyeron covariables estacionales en el ajuste, la fecha de fin es completamente arbitraria. Caso contrario, la fecha de comienzo no puede ser anterior al inicio de la serie de covariables, tampoco puede ser posterior. Se debe introducir una fecha en formato datecontrol: objeto de control creado con la funcióncontrol_sim().output_folder: ruta al directorio donde se guardarán los resultados, tanto finales como intermedios.output_filename: nombre del archivo de salida. Para facilitar la interoperabilidad, el archivo generado es un archivo de texto en formato separado por comas (.csv)seasonal_covariates: datos agregados trimestrales. Si el ajuste se realizó con covariables, la generación también debe realizarse con ellas. Caso contrario se producirá un error. Se debe introducir un data frame con los valores agregados para las tres variables (precipitación y temperaturas máxima y mínima) pero no necesariamente deben ser los mismos a los utilizados en el ajuste. Si se desean simular tendencias de algún tipo, ya sea de un modelo de cambio climático o arbitrarias, se deben perturbar estas variables trimestrales e introducirlas aquí. Estas series si deben tener la misma longitud que el período a generarverbose: controla la impresión de mensajes en la consola. FALSE por defecto.
# Al correr la función se realiza la generación de series para cada una de las estaciones.
# En este caso, por cuestiones de tiempo, vamos a cargar un objeto con los resultados de la simulación
simulated_climate <- gamwgen::spatial_simulation(model = gamgen_fit, # Objeto con los resultados del ajuste
simulation_locations = grilla_simulacion_centers, # Estaciones para las cuales simular
start_date = as.Date('2019-01-01'), # Fecha de comienzo de las simulaciones
end_date = as.Date('2019-01-31'), # Fecha de fin de las simulaciones
control = control_sim, # Objeto con la configuración
output_folder = getwd(), # Directorio donde se guardarán los resultados
output_filename = 'simulations.csv', # Nombre del archivo de salida
seasonal_covariates = seasonal_covariates, # Covariables estacionales
verbose = FALSE) # Impresión de mensajes en la consolaEsta función produce dos tipos de resultados: una lista que permanece en el ambiente de R y los datos generados que son guardados como .csv en el directorio indicado precedentemente.
# Copiamos el archivo preajustado a nuestro directorio de trabajo
if (!fs::file_exists('output_data/spatial/simulated_spatial_conditional.RData')) {
fs::file_copy(system.file('/autorun/spatial', "simulated_spatial_conditional.RData", package = "gamwgen"),
new_path = 'output_data/spatial/simulated_spatial_conditional.RData')
}
# Cargamos el archivo recientemente creado
load('output_data/spatial/simulated_spatial_conditional.RData')
# Clase del objeto con el ajuste del generador
class(simulated_climate)## [1] "list" "gamwgen.climate"## [1] "nsim"
## [2] "seed"
## [3] "realizations_seeds"
## [4] "simulation_points"
## [5] "output_file_with_results"
## [6] "output_file_fomart"
## [7] "rdata_file_with_fitted_stations_and_climate"
## [8] "exec_times"La lista contiene los siguientes objetos:
nsim: cantidad de realizaciones.seed: semilla general para toda la generación. Corresponde a la que se incluye en la función de control.realization_seeds: semillas para cada una de las realizaciones. Esto permite replicar los resultados.simulation_points: puntos donde se generaron las series sintéticas.output_fil_with_results: nombre del archivo con los resultados.output_file_format: tipo de archivo de salida, en este caso .csv.rdata_file_with_fitted_stations_and_climate: archivo.RDatacon los datos meteorológicos observados que fueron utilizados en el ajuste. También se incluyen los metadatos de cada uno de esos puntos.exec_times: tiempo de ejecución del la generación.
Ahora veremos el formato del archivo de salida que contiene las series sintéticas.
Para desmenuzar la generación del tiempo local se pueden correr algunas de las funciones que forman parte de local_simulation. Por ejemplo, del objeto gamgen_fit se extraen los residuos y con la función gamwgen:::generate_residuals_statistics se calculan los parámetros.
En la tabla anterior se muestran los parámetros necesarios para generar el tiempo local de las temperaturas máxima y mínima.
gen_noise_params <- gamwgen:::generate_month_params(
residuals = gamgen_fit$models_residuals,
observed_climate = gamgen_fit$models_data,
stations = gamgen_fit$stations)
# Ejemplo de variograma de residuos de temperatura máxima para días secos
gen_noise_params[[1]]$variogram_parameters$tmax_dry ## [1] 0.000000 7.045242 493.298928Explicar variograma.
Los parámetros para cada uno de los meses permiten generar valores de tiempo local para las dos temperaturas a partir de una distribución normal multivariada.
A continuación se muestra un ejemplo para el año 2019 sólo para los días lluviosos.
RandomFields::RFoptions(printlevel = 0, spConform = FALSE)
grilla_simulacion_centers <- uruguay %>%
sf::st_transform(gamgen_fit$crs_used_to_fit) %>%
sf::st_make_grid(cellsize = 25000, what = 'centers') %>%
sf::st_sf(grid_id = 1:length(.), geometry = .) %>%
sf::st_intersection(uruguay) %>%
dplyr::select(geometry) %>%
dplyr::mutate(latitude = sf::st_coordinates(.)[,2],
longitude = sf::st_coordinates(.)[,1])
# Mapa de Uruguay
uruguay <- raster::getData('GADM', country='URY', level=1) %>%
sf::st_as_sf(.) %>%
sf::st_transform(crs = 32721)
temp_local <- control_sim$temperature_noise_generating_function(
simulation_points = grilla_simulacion_centers,
gen_noise_params = gen_noise_params,
month_number = 1L,
selector = 'Dry',
seed = 1234) %>%
dplyr::mutate(date = as.Date(paste0('2019-', '01', '-', '01'))) %>%
dplyr::filter(date == as.Date('2019-01-01')) %>%
tidyr::gather(variable, valor, -c('geometry'))
tmax_residuals <- gamwgen:::sf2raster(temp_local %>%
dplyr::filter(variable == 'tmax_residuals'),
'valor') %>%
raster::as.data.frame(., xy = TRUE) %>%
dplyr::rename(., 'longitude' = 'x', 'latitude' = 'y') %>%
dplyr::mutate(., x = longitude, y = latitude) %>%
sf::st_as_sf(., coords = c( 'x', 'y')) %>%
sf::st_set_crs(., 32721) %>%
sf::st_intersection(uruguay) %>%
dplyr::select(latitude, longitude, z, geometry)
colr <- colorRampPalette(RColorBrewer::brewer.pal(9, 'YlOrRd'))
quiebres <- c(0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0)
etiquetas <- c('0', '5', '10', '15', '20', '25', '30', '35', '40', '45')
ggplot(data = tmax_residuals,
ggplot2::aes(x = longitude, y = latitude, fill = z)) +
geom_tile(colour="black") +
scale_fill_gradientn(name = "Temperatura máxima [°C]",
colours = colr(9),
breaks = quiebres,
labels = etiquetas) +
ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE) +
ggplot2::theme_bw() +
ggplot2:: labs(title = "", x = "", y = "") +
ggplot2::theme(panel.background = element_rect(fill = "ivory"),
axis.text = element_text(size = 11, colour = 1))
La líonea naranja corresponde a la serie para temperaturas máximas y la azul para temperaturas mínimas.
# Se carga el set de datos simulados
simulated_climate <- readr::read_csv(here::here('output_data/spatial/simulated_spatial_conditional.csv')) ##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## realization = col_double(),
## point_id = col_double(),
## longitude = col_double(),
## latitude = col_double(),
## date = col_date(format = ""),
## tmax = col_double(),
## tmin = col_double(),
## prcp = col_double()
## )| realization | point_id | longitude | latitude | date | tmax | tmin | prcp |
|---|---|---|---|---|---|---|---|
| 1 | 1 | 453759.4 | 6590418 | 2019-01-01 | 33.03589 | 19.72844 | 1.826571 |
| 1 | 2 | 553759.4 | 6590418 | 2019-01-01 | 31.77257 | 19.65962 | 27.883864 |
| 1 | 3 | 578759.4 | 6590418 | 2019-01-01 | 32.71455 | 20.47479 | 32.879767 |
| 1 | 4 | 428759.4 | 6615418 | 2019-01-01 | 30.64363 | 16.02136 | 7.032041 |
| 1 | 5 | 453759.4 | 6615418 | 2019-01-01 | 31.83107 | 18.12258 | 14.671730 |
| 1 | 6 | 478759.4 | 6615418 | 2019-01-01 | 33.32291 | 18.75998 | 14.781878 |
| 1 | 7 | 503759.4 | 6615418 | 2019-01-01 | 32.22893 | 18.78158 | 11.007127 |
| 1 | 8 | 528759.4 | 6615418 | 2019-01-01 | 31.95957 | 17.71497 | 24.732878 |
| 1 | 9 | 553759.4 | 6615418 | 2019-01-01 | 31.93247 | 17.65988 | 26.055372 |
| 1 | 10 | 453759.4 | 6640418 | 2019-01-01 | 30.67991 | 17.13227 | 16.806296 |
El resultado de la generación es un archivo .csv que contiene la siguiente información:
realization: número de realización. Es un valor entero entre 1 y la cantidad de realizaciones definida por el usuario.station_id: número unívoco de identificación de la estación meteorológica o del punto arbitrario.date: fechas de cada uno de los días de la simulación.tmax: valores de temperatura máxima generada expresada en °C.tmin: valores de temperatura mínima generada expresada en °C.prcp: valores de precipitación diaria generada expresada en mm.
La siguiente Figura muestra un ejemplo de las series de temperaturas máximas y mínimas generadas.
uruguay <- raster::getData('GADM', country='URY', level=1) %>%
sf::st_as_sf(.) %>%
sf::st_transform(crs = 32721)
simulated_climate_geo <- simulated_climate %>%
dplyr::filter(date == as.Date('2019-01-01'),
realization == 1) %>%
sf::st_as_sf(coords = c('longitude', 'latitude'), crs = 32721) %>%
gamwgen:::sf2raster(., 'prcp') %>%
raster::as.data.frame(., xy = TRUE) %>%
dplyr::rename(., 'longitude' = 'x', 'latitude' = 'y') %>%
dplyr::mutate(., x = longitude, y = latitude) %>%
sf::st_as_sf(., coords = c( 'x', 'y')) %>%
sf::st_set_crs(., 32721) %>%
sf::st_intersection(uruguay) %>%
dplyr::select(latitude, longitude, z, geometry) ## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriescolores <- colorRampPalette(c("#cccccc", "#7bccc4", "#2b8cbe", "#0868ac",
"#084081", "#807dba", "#6a51a3", "#54278f", "#3f007d", "#e7298a",
"#ce1256", "#67001f"))
quiebres <- c(0.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0)
etiquetas <- c('0', '5', '10', '20', '30', '40', '50', '60', '70', '80', '90', '100', '110', '120', '130', '140', '150', '160')
ggplot2::ggplot(data = simulated_climate_geo) +
ggplot2::geom_tile(ggplot2::aes(x = longitude, y = latitude, fill = z)) +
ggplot2::scale_fill_gradientn(name = "Precipitacópn [mm]",
colours = colores(17),
breaks = quiebres,
labels = etiquetas,
limits = c(0, 100)) +
ggplot2::geom_sf(data = uruguay, fill = NA, color = "black", inherit.aes = FALSE)
#cowplot::ggdraw() +
# cowplot::draw_plot(tmax_plot, x = 0, y = 0.5, width = 1, height = .5) +
# cowplot::draw_plot(tmin_plot, x = 0, y = 0, width = 1, height = .5)En el panel superior se muestra la temperatura máxima observada (negra) y las temperaturas sintéticas (naranja) para el año 2019. En el infeior se muestra la temperatura mínima diaria observada (negra) y cada una de las realizaciones (azul).
Agregar covariables interpoladas!!!!
Modelo no lineal
Modelos Lineales (Generalizados)
GLM: terminos lineales GAM: terminos no necesariamente lineales
Considerando \(n\) observaciones, \(x_{i}\), \(y_{i}\), donde \(y_{i}\) es una observación de una variable aleatoria, \(Y_{i}\), con una esperanza de \(\mu_{i}\) ≡ \(E(Y)\). Un modelo apropiado para expresar la relación entre \(x\) e \(y\) es el siguiente:
\[ y_i = \mu_i + \epsilon_i \] donde, \(\mu_i\) = \(x_{i}\beta\)
\(\beta\), parámetro desconocido, y \(\epsilon_i\) son variables aleatorias con media cero y cada una con variancia \(\sigma^2\). Por lo que el modelo establece que \(Y\) está dado por \(x\) multiplicada por una constante más un termino aleatorio. \(Y\) es un ejemplode variable respuesta mientras que \(x\) es un ejemplo de un predictor.
El modelo lineal se describe mediante la siguiente ecuación:
\[ y_i = \beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2 + \ldots + \epsilon_i \] \(y_i\) corresponde a la respuesta como una combinación lineal de covariables, \(x_{ji}\), con un offset.
La variable \(y_i\sim\) puede ser de cualquier distribución de la familia exponencial (Normal, Poisson, etc.)
El termino de error \(\epsilon_i\) está linealmente distribuido
Una vez obtenida la cantidad de simulaciones deseadas, se procede a la validación, que consiste en contrastar dichos valores con los registros observados para evaluar cuán bien éstas han capturado las propiedades estadísticas de los registros históricos. Por último, una vez que los resultados son validados, es posible simular series sintéticas con confianza tanto en localidades, con o sin datos históricos, como sobre una grilla regular que represente a una región.
El modelo es capaz de generar dos tipos de series sintéticas: no condicionadas y condicionadas. Las primeras, son series estacionarias, es decir, que no tienen tendencia temporal por lo que representan las secuencias medias de una localidad o región. Las series condicionadas, en cambio, incluyen en su proceso de ajuste variables auxiliares que reproducen cambios observados (ej., variabilidad climática natural) o proyectados en el clima. Estas variables auxiliares (también llamadas “covariables”) son totales trimestrales de lluvia acumulada y medias trimestrales de temperatura máxima y mínima. Su inclusión permite modelar cambios interdecadales o probabilidades estacionales, por lo que se pueden utilizar salidas de pronósticos estacionales, de modelos de cambio climático, o incluso valores arbitrarios definidos por el usuario para crear escenarios de interés (por ejemplo, un aumento sostenido de la precipitación). Si en cambio se utilizan covariables calculadas a partir de los datos, se obtienen series que hemos denominado pseudohistóricas.
Este es un ejemplo de aplicación del Generador Estocástico diario y multisitio basado en GAMs. Para ello se ajustará el modelo local del generador para una localidad de la Pampa Húmeda Argentina.
Is this linear?
set.seed(2) ## simulate some data...
dat <- gamSim(1, n = 400, dist = "normal", scale = 0.5, verbose=FALSE)
dat <- dat[,c("y", "x0", "x1", "x2", "x3")]
p <- ggplot(dat,aes(y=y,x=x1)) +
geom_point() +
theme_minimal()
print(p)
References
Kleiber, W., Katz, R.W., Rajagopalan, B., 2012. Daily spatiotemporal precipitation simulation using latent and transformed gaussian processes. Water Resources Research 48.
Kleiber, W., Katz, R.W., Rajagopalan, B., 2013. Daily minimum and maximum temperature simulation over complex terrain. Ann. Appl. Stat. 7, 588–612.
Simpson, G., 2018. Introducing gratia. From the Bottom of the Heap.
Verdin, A., 2016. Stochastic space-time modeling for agricultural decision support in the argentine pampas (Thesis).
Wickham, H., 2011. Ggplot2. Wiley Interdisciplinary Reviews: Computational Statistics 3, 180–185.
Wood, S.N., 2001. Mgcv: GAMs and generalized ridge regression for r. R news 1, 20–25.
Wood, S.N., 2017. Generalized additive models: An introduction with r. Chapman; Hall/CRC.